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Abstract 

Dual  photography  is  a  technique  capable  of  reconstructing  an  image  of  an  occluded 
scene  from  reflected  light  by  exploiting  Helmholtz  reciprocity.  The  primary  limitation 
of  dual  photography  is  the  line-of-sight  requirement  of  the  source  used  to  illuminate 
the  occluded  scene.  Complex  radiometric  modeling  of  multiple  reflections  allowed  a 
technique  called  indirect  photography  to  overcome  the  line-of-sight  requirement  of 
dual  photography  and  recover  some  information  from  the  hidden  scene  that  was  not 
directly  visible  from  the  camera  or  the  light  source. 

This  research  focuses  on  reflective  inverse  diffusion,  which  was  a  proof-of-concept 
experiment  that  used  phase  modulation  to  shape  the  wavefront  of  a  laser  causing  it 
to  refocus  after  reflection  from  a  rough  surface.  By  refocusing  the  light,  reflective 
inverse  diffusion  has  the  potential  to  eliminate  the  complex  radiometric  model  of 
indirect  photography  by  creating  a  virtual  light  source  at  the  first  diffuse  reflector 
that  satisfies  the  line-of-sight  requirement  of  dual  photography.  However,  the  initial 
reflective  inverse  diffusion  experiments  provided  no  mathematical  background  and 
were  conducted  under  the  premise  that  the  process  operated  similarly  to  transmissive 
inverse  diffusion. 

In  this  research,  diffraction  modeling  of  the  reflective  inverse  diffusion  experiments 
led  to  the  development  of  Fourier  transform-based  simulations.  Simulations  and  ex¬ 
perimentation  were  used  to  develop  reflection  matrix  methods  that  determine  the 
proper  phase  modulation  to  refocus  light  after  reflection  to  any  location  in  the  ob¬ 
servation  plane.  These  techniques  provide  a  new  method  for  controlled  illumination 
of  an  occluded  scene  that  can  be  used  in  conjunction  with  dual  photography.  This 
document  provides  the  mathematical  background  for  reflective  inverse  diffusion,  the 
reflection  matrix  methods  for  phase  modulation,  and  describes  the  simulations  and 
experiments  conducted. 
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REFLECTION  MATRIX  METHOD  FOR  CONTROLLING  LIGHT  AFTER 
REFLECTION  FROM  A  DIFFUSE  SCATTERING  SURFACE 

I.  Introduction 

The  Helmholtz  Reciprocity  Theorem  [1]  describes  the  symmetric  nature  of  light 
propagation.  This  theorem  states  that  the  effect  measured  at  location  P  from  a 
point  source  at  location  P0  would  be  the  same  as  the  effect  measured  at  P0  from 
an  equivalent  point  source  located  at  P.  Dual  Photography  [2]  is  a  technique  that 
takes  advantage  of  Helmholtz  reciprocity.  The  process  involved  a  scene  illuminated 
by  a  digital  projector  pixel  by  pixel.  A  digital  camera  then  recorded  the  reflected 
light  from  the  scene  for  each  pixel  illumination.  Using  this  information,  a  transport 
matrix  was  produced  to  relate  the  effects  of  the  projector  pixels  to  the  camera  pixels. 
The  transport  matrix  can  be  used  to  mathematically  interchange  the  positions  of 
the  camera  and  the  projector  producing  a  view  of  the  scene  from  the  position  of  the 
projector  as  if  it  were  illuminated  from  the  position  of  the  camera.  This  allowed 
recovery  of  scene  information  that  was  only  visible  from  the  projector’s  point  of 
view  [2], 

Scene  information  hidden  from  the  camera  but  recovered  using  dual  photography 
must  be  visible  from  the  viewpoint  of  the  projector.  This  requirement  limits  its 
usefulness  outside  the  author’s  original  work  of  capturing  and  relighting  scenes.  In 
other  words,  if  the  desired  image  is  observable  from  the  viewpoint  of  the  projector, 
then  the  simplest  solution  would  be  to  place  a  camera  there.  Indirect  photography  [3, 
4]  is  a  method,  developed  at  the  Air  Force  Institute  of  Technology  (AFIT),  to  remove 
the  line-of-sight  requirement  of  dual  photography  and  recover  scene  information  not 
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directly  visible  from  the  position  of  the  camera  or  projector.  This  would  allow  the  light 
source,  now  a  laser,  and  the  camera  to  be  co-located  and  yet  still  retain  the  ability 
“to  see  around  corners.”  This  was  accomplished  through  radiometric  modeling  of 
multiple  diffuse  reflections.  Indirect  photography  was  successful  at  recovering  some 
hidden  scene  information  after  two  diffuse  reflections.  Initial  methods  for  indirect 
photography  limited  reflections  to  forward  scattering  only  [3,4], 

Indirect  photography  illuminates  a  scene  with  reflectively  scattered  light.  Since 
the  illumination  pattern  no  longer  provides  a  controlled  canonical  basis,  as  with  dual 
photography,  the  ability  to  construct  a  proper  transport  matrix  to  image  the  occluded 
scene  is  significantly  impaired.  Adaptive  optics  techniques,  such  as  phase  modulation, 
can  be  implemented  to  compensate  for  the  reflective  scattering  of  the  illumination 
source.  The  ability  to  control  the  illumination  source  and  refocus  light  after  reflection 
would  reduce  the  indirect  photography  problem  back  to  that  of  dual  photography  and 
allow  proper  construction  of  a  transport  matrix  for  imaging  the  hidden  scene. 

Transmissive  inverse  diffusion  [5-8]  is  a  method  for  refocusing  light  after  being 
scattered  by  transmission  through  a  turbid  medium.  Reflective  inverse  diffusion  is 
developed  in  this  research.  The  details  of  reflective  inverse  diffusion  will  be  discussed 
in  Chapter  II.  Applying  the  concept  of  inverse  diffusion  to  indirect  photography, 
refocusing  light  from  the  first  diffuse  reflection,  has  the  potential  to  simplify  the 
method  back  to  that  of  dual  photography.  From  the  observation  plane  perspective, 
the  origin  of  the  light  source  is  the  first  diffuse  reflector,  thus  satisfying  the  line-of- 
sight  requirement  of  dual  photography. 

Reflective  inverse  diffusion  was  a  proof-of-concept  research  at  AFIT.  Phase  mod¬ 
ulation  was  employed  to  shape  a  plane  wave  prior  to  reflecting  off  a  diffuse  surface. 
Phase  modulation  was  achieved  via  a  liquid  crystal  on  silicon  (LCoS)  spatial  light 
modulator  (SLM).  The  results  showed  that  a  properly  modulated  plane  wave  could 
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be  focused  to  a  tight  spot  after  reflection.  The  SLM  also  demonstrated  some  degree 
of  control  over  the  location  of  the  focus  spot  in  the  observation  plane  [9]. 

The  purpose  of  this  research  is  to  develop  a  method  of  controlled  illumination  using 
reflectively  scattering  light  that  could  be  used  in  conjunction  with  dual  photography 
to  allow  imaging  around  corners.  This  method  required  significant  advancement  of 
the  theory  and  application  of  reflective  inverse  diffusion.  The  proof-of-concept  ex¬ 
periments  were  based  on  methods  developed  for  transmissive  inverse  diffusion.  Both 
transmissive  inverse  diffusion  and  the  original  experimental  results  for  reflective  in¬ 
verse  diffusion  are  covered  in  Chapter  II.  Diffraction-based  models  were  developed 
to  provide  the  basis  for  simulations  of  reflection  inverse  diffusion  and  are  presented 
in  Chapter  III.  The  reflection  matrix  (RM)  was  developed  to  provide  the  necessary 
control  of  the  reflected  light.  The  relationship  between  every  segment  of  the  SLM 
and  every  segment  of  the  charge-coupled  device  (CCD)  detector  is  contained  in  the 
RM.  The  algorithm  used  to  measure  and  construct  the  RM  is  covered  in  Chapter  IV. 
The  experimental  optical  setup  and  reflector  material  surface  properties  that  effect 
the  performance  of  the  RM  are  examined  in  Chapter  V.  Finally,  initial  methods  of 
reconstructing  an  RM  from  partial  data  and  expanding  an  existing  RM  to  cover  an 
observation  plane  larger  than  that  measured  are  presented  in  Chapter  VI. 
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II.  Background 


The  concept  of  reflective  inverse  diffusion  was  based  on  methods  developed  for 
refocusing  light  through  thin  films  of  a  turbid  medium.  When  coherent  light  is  trans¬ 
mitted  through  a  stationary  diffuser  (i.e.  a  turbid  medium),  a  fine  granular  intensity 
pattern  is  formed  called  speckle  [10].  The  speckle  is  caused  by  the  individual  path 
length  differences  imparted  by  the  medium  to  different  portions  of  the  coherent  light. 
Thus,  transmissive  inverse  diffusion  uses  phase  modulation  to  conjugate  the  phase 
changes  imparted  by  the  turbid  medium  and  causes  the  light  to  refocus  after  trans¬ 
mission. 

Speckle  can  also  be  produced  from  coherent  light  reflecting  off  a  rough  surface. 
The  path  length  differences  that  produce  the  speckle  are  caused  by  the  surface  height 
deviations  of  the  material  [10].  A  rough  surface,  in  the  context  of  optics,  is  any 
material  with  a  surface  height  standard  deviation  on  the  order  of  a  single  wavelength 
of  light  (A)  [11].  Due  to  the  similarities  between  transmissive  and  reflective  speckle, 
reflective  inverse  diffusion  methods  were  initially  based  on  the  phase  modulation 
techniques  used  in  transmissive  inverse  diffusion. 

2.1  Transmissive  Inverse  Diffusion 

Driven  by  potential  medical  applications  of  imaging  through  tissue,  inverse  diffu¬ 
sion  was  first  accomplished  in  transmission  by  focusing  light  through  a  turbid  medium. 
Vellekoops’s  method  exploited  the  linearity  of  the  scattering  process  [5].  The  exper¬ 
imental  set-up  for  transmissive  inverse  diffusion  is  shown  in  Figure  1.  The  turbid 
medium  is  placed  between  two  microscope  objectives.  The  first  microscope  objective 
applies  a  demagnified  image  of  the  modulated  beam  from  the  SLM  onto  the  turbid 
medium.  The  second  objective  collects  the  transmitted  scatter  pattern  and  provides 
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a  magnified  image  on  the  CCD  detector  that  provides  feedback  for  the  system. 

Both  the  SLM  and  CCD  are  addressed  in  segments  and  each  segment  consists  of 
multiple  pixels.  The  basic  iterative  algorithm  for  transmissive  inverse  diffusion  selects 
one  segment  of  the  CCD  and  seeks  to  maximize  the  intensity  at  that  location.  The 
phase  of  each  of  the  SLM  segments  is  cycled  through  a  small  subset  of  the  available 
[0,  27 r]  range  until  the  maximum  intensity  at  the  CCD  target  segment  is  achieved. 


Microscope 

Objectives 


Figure  1.  Optical  setup  for  transmissive  inverse  diffusion.  A  microscope  objective  is 
used  to  demagnify  the  modulated  beam  and  image  it  onto  the  scattering  medium.  The 
transmitted  scatter  pattern  is  magnified  by  a  second  microscope  objective  and  image 
onto  the  feedback  CCD. 


The  transmitted  field  at  the  mth  segment  of  the  observation  plane  is  given  by, 

N 

F  —  f  A  p^n  Mj 

J-Jrri  /  l,mnnnc  V  ^  / 

n= 1 

where  tmn  represents  the  mnth  element  of  the  complex-valued  transmission  matrix, 
and  Ane^n  is  the  field  from  the  nth  segment  of  the  SLM.  The  source  field  is  segmented 
by  the  modulator  into  N  segments  and  transmitted  through  the  scattering  medium. 
The  field  in  the  observation  plane  consists  of  a  linear  combination  of  these  N  segments 
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of  the  source  field  at  each  of  the  m  observation  positions  [5].  Normalizing  equation 
(1)  to  intensity,  An  =  1/Vn)  and  the  observed  intensity  at  the  mth  segment  is  then 


n=  1 


(2) 


Intensity  enhancement  is  a  measure  of  performance  for  transmissive  inverse  diffu¬ 
sion  and  was  defined  by  Vellekoop  as  [7], 


n  = 


{I opt) 


(3) 


\Jrnd) 

which  is  a  simple  ratio  of  the  average  intensity  of  the  optimized  spot  {I opt)  divided 
by  the  average  intensity  of  the  unoptimized  random  speckle  (Irnd)  [6].  Each  grain  of 
speckle  in  the  intensity  pattern  is  the  sum  of  a  large  number  of  individual  light  paths 
through  the  medium.  Applying  the  Central  Limit  Theorem,  the  distribution  of  the 
individual  paths  through  the  medium  is  approximated  by  a  Gaussian  [12].  Assuming 
the  paths  through  the  material  are  independent,  the  real  and  imaginary  components 
of  the  turbid  medium  become  individually  Gaussian,  then  this  distribution  is  known 
to  be  a  circular  Gaussian  [13].  Thus,  for  transmission  through  a  random  medium,  the 
tmn  terms  follow  a  statistically  independent  complex  circular  Gaussian  distribution 
[13]  with  properties  that  can  be  used  to  simplify  equation  (3)  and  express  the  ideal 
intensity  enhancement,  r)transmissive,  as  a  function  of  the  total  number  of  SLM  segments 

[7], 


7 T  7 r 

Vtranmissive  ^  (  A  1)  T  1  ~  A7 


(4) 
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2.2  Transmission  Matrices 


Matrix  multiplication  has  been  used  to  represent  the  linear  combination  of  N 
segments  of  the  modulator  from  Vellekoop’s  method.  The  transmission  matrix.  T, 
approximates  the  effects  of  scattering  by  the  medium  and  propagation  through  the 
system.  Several  methods  for  experimentally  determining  the  transmission  matrix  all 
rely  on  orthogonal  Hadamard  or  Walsh- Hadamard  basis  vectors  [14-16].  Recently, 
methods  for  determining  transmission  matrices  using  the  temporal  Fourier  transfrom 
of  instensity  measurements  have  been  developed  [17].  The  advantage  of  the  matrix 
approach  is  the  ability  to  use  the  transmission  matrix  to  determine  the  input  field  for 
any  desired  output. 

2.2.1  Other  Literature. 

Iterative  methods  for  transmissive  inverse  diffusion  were  developed  by  Vellekoop 
and  Mosk  [7].  In  general,  these  methods  show  a  linear  improvement  in  the  enhance¬ 
ment  of  the  tightly  focused  spot  with  an  increase  in  the  number  of  segments  of  the 
modulator.  Optimization  for  a  single  SLM  segment  took  1.2  seconds  [5].  However, 
it  is  estimated  that  a  single  segment  optimization  of  1  millisecond  is  needed  to  make 
dynamic  measurements  in  biological  tissue  [5-8,18].  This  optimization  time  has  been 
achieved  by  using  acousto-optic  modulators  to  rapidly  scan  for  the  optimum  phase 
before  applying  it  to  the  relatively  slow  SLM.  Feedback  time  was  reduced  by  us¬ 
ing  a  simple  photo-diode  for  intensity  measurements  in  the  observation  plane  [19]. 
Phase  modulation  has  also  been  used  to  maximize  the  transmission  of  light  through 
scattering  media  without  focusing  it  to  a  single  spot  [20] . 

Transmission-matrix  approaches  were  soon  developed  beyond  the  iterative  meth¬ 
ods  [15,21],  These  methods  tend  to  be  favored  over  iterative  methods  as  they  allow 
for  control  of  the  light  in  the  observation  plane  without  need  to  rerun  a  sometimes 
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lengthy  iterative  algorithm.  Several  methods  for  measuring  the  transmission  matrix 
involve  using  the  Hadamard  orthonormal  basis  vectors  to  interrogate  the  scattering 
sample  [15,16,21];  however,  more  recent  publications  allow  for  any  orthogonal  basis 
to  be  used  [17].  Statistics  and  control  capabilities  of  transmission  matrices  have  been 
examined  in  detail  [22],  Most  wavefront-shaping  methods  favor  phase  modulation; 
however,  methods  for  determining  the  transmission  matrix  of  a  scattering  medium 
have  been  developed  using  binary  amplitude  modulation  with  digital  micro- mirror  de¬ 
vices  (DMDs).  These  devices  operate  in  the  20  kHz  range  much  faster  than  the  30-60 
Hz  of  liquid  crystal  on  silicon  (LCoS)  spatial  light  modulators  (SLMs)  [14,23].  Back¬ 
ground  speckle  tends  to  be  higher  for  images  generated  using  transmission  matrices 
due  to  errors  in  transmission  matrix  measurement  and  estimation  process;  however, 
there  are  published  methods  for  reducing  this  side  effect  [24] . 

Transmissive  inverse  diffusion  has  been  used  to  image  through  scattering  ma¬ 
terial  [25,26].  Imaging  through  biological  tissues  along  with  other  applications  in 
biophotonics  have  also  been  explored  [27].  Transmission  matrix  methods  have  also 
been  used  to  determine  the  polarization  state  of  incident  light  [28] .  Wavefront  shaping 
methods  using  phase  modulation  were  used  to  create  a  programmable  optical  circuit 
using  scattering  materials  [29].  The  ability  to  control  the  propagation  time  [30]  and 
change  the  direction  of  propagation  in  scattering  media  [31]  has  also  been  demon¬ 
strated.  Both  Mosk  and  Vellekoop  recently  published  summaries  covering  most  of 
the  new  techniques  for  transmissive  inverse  diffusion  using  wavefront  shaping  [32,33], 
where  it  was  mentioned  that  the  transmission  matrix  approach  could  be  applied  to 
any  linear  process  without  loss  of  generality,  including  reflection. 


2.3  Reflective  Inverse  Diffusion 


2.3.1  Experimental  Set-up. 

The  experimental  set-up  for  reflective  inverse  diffusion  was  adapted  from  that  of 
transmissive  inverse  diffusion  shown  in  Figure  1.  Although  similar,  the  two  have 
significant  differences.  In  transmission,  microscope  objectives  are  located  on  both 
sides  of  the  scattering  medium.  To  have  optical  elements  this  close  to  the  scattering 
medium  is  not  practical  for  use  in  reflection,  and  is  also  an  impractical  requirement 
for  imaging  around  corners  [9]. 


SLM 


Figure  2.  Focal  plane  optical  setup  for  reflective  inverse  diffusion.  A  vertically  polarized 
HeNe  laser  is  expanded,  collimated,  and  normally  incident  on  the  SLM.  The  phase 
modulated  beam  is  then  focused  onto  the  rough  surface  reflector  with  positive  lens 
(L i)  and  the  reflected  speckle  pattern  is  recorded  by  the  CCD.  The  mirror  (Mi)  and 
the  non-polarizing  beam  splitter  (NPBS)  are  used  to  allow  normal  incidence  with  the 
SLM.  For  focal-plane  experiments  and  simulations,  the  lens  focal  length,  /,  is  500  mm, 
and  the  distances  Z\  and  Z2  are  15  ±  0.5  cm  and  50  ±  0.5  cm,  respectively. 


The  original  proof-of-concept  experiment  for  reflective  inverse  diffusion  (shown  in 
Figure  2)  was  done  with  a  5  mW  HeNe  laser  (632.8  nm  wavelength)  linearly  polarized 
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with  the  vertical  axis  of  the  SLM.  The  beam  was  then  expanded  and  collimated.  A 
non-polarizing  beam  splitter  (NPBS)  allows  normal  incidence  onto  the  SLM.  The 
modulated  beam  was  focused  onto  the  diffuse  reflector  using  a  500-mm  positive  lens. 
The  focused  beam  is  incident  on  the  reflector  at  45°,  the  reflector  then  scatters  the 
beam  onto  the  CCD  detector  [9]. 

2.3.2  Algorithm. 

The  algorithm  for  reflective  inverse  diffusion  was  adapted  from  transmissive  inverse 
diffusion,  a  method  for  focusing  light  after  transmission  through  a  diffuse  media.  The 
basic  algorithm  measures  the  intensity  in  the  observation  plane  while  incrementally 
adjusting  the  phase  of  the  SLM  pixel.  The  phase  value  that  produced  the  highest 
intensity  spot  in  the  observation  plane  was  selected,  and  the  algorithm  moved  to  the 
next  pixel  [5,9]. 

The  SLM,  used  in  the  experiment,  was  a  Boulder  Nonlinear  Systems  (BNS)  model 
P512  with  a  resolution  of  512-by-512  pixels.  Each  pixel  was  15/irn  square,  and  capable 
of  256  discrete  phase  levels  with  a  total  phase  stroke  of  2n.  The  SLM  is  therefore 
capable  of  226  (i=s  108)  different  phase  screens.  Since  it  is  impractical  to  test  all 
possible  pixel/phase  combinations,  the  SLM  was  further  grouped  into  16-by-16  pixels 
blocks  to  be  modulated  as  a  single  super  pixel;  the  resolution  of  SLM  would  then 
be  32  by  32  super  pixels.  Each  super  pixel  was  phase  modulated  from  0  to  2n  in 

7T 

increments  of  —  for  a  total  of  21  possible  phase  values  [9]. 

The  algorithm  is  forward  only;  once  a  super  pixel  is  optimized,  it  moves  to  the 
next  super  pixel.  The  algorithm  never  determines  if  the  new  phase  value  for  the 
current  super  pixel  affects  any  of  the  previously  optimized  super  pixels.  Since  the 
SLM  is  phase  only,  the  effects  of  modulating  a  single  super  pixel  are  not  independent 
of  each  other.  It  can  be  shown  that  with  a  single  pass  through,  the  algorithm  does 
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not  produce  a  unique  solution,  and  that  the  solution  depends  on  the  initial  values  of 
the  SLM  pixels. 

Looping  through  the  algorithm  does  improve  spot  intensity,  and  reduces  back¬ 
ground  speckle;  however,  it  is  very  time  consuming.  Coarse  pre-optimization  was 
employed  as  a  compromise  between  spot  enhancement  and  processing  time.  The 
SLM  is  initially  segmented  into  2-by-2  blocks,  each  consisting  of  256-by-256  pixels. 
Each  block  is  optimized  to  one  of  the  21  phase  values.  The  SLM  is  then  be  seg¬ 
mented  into  4-by-4  blocks  and  optimized  to  one  of  the  21  phase  values.  This  process 
of  course  optimization  is  repeated  until  the  SLM  is  segmented  into  the  32-by-32  super 
pixels  [5,9]. 

The  initial  algorithm  sought  to  improve  efficiency  by  under-filling  the  SLM,  and 
only  modulating  the  pixels  illuminated  by  the  laser.  This  only  produced  marginal 
reductions  in  processing  time,  but  added  the  uncertainty  of  laser  alignment.  It  would 
be  impossible  to  determine  if  some  SLM  pixels  were  illuminated  by  the  laser,  but  left 
unmodulated.  The  end  result  of  the  algorithm  was  1,160  super  pixel  optimizations 
for  each  experiment.  Each  super-pixel  optimization  takes  approximately  32  seconds 
to  process.  The  total  process  took  over  10  hours  to  complete  [9]. 

2.3.3  Summary  of  Results. 

The  original  proof-of-concept  experiments  for  reflective  inverse  diffusion  measured 
the  enhancement  from  six  different  scattering  materials:  Spectralon®,  brushed  alu¬ 
minum,  sandblasted  aluminum,  Infragold®,  white  paint  on  glass,  and  graphite.  The 
reflector  materials  were  selected  based  on  the  differences  in  scattering  properties.  The 
surface  roughness  of  the  samples  was  measured  and  compared  with  the  reflected  spot 
enhancement,  and  the  final  full-width-at-half-maximum  (FWHM)  spot  size  produced 
by  the  algorithm.  These  measurements  are  shown  in  Table  1. 
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Table  1.  Summary  values  of  enhancement,  surface  roughness,  and  spot  size  for  each  of 
the  six  reflective  samples  measured  and  the  one  transmission  sample. 


Reflective  Samples 

Enhancement  (rj) 

Roughness  ( Rrms ) 

Final  FWHM 

Brushed  Aluminum 

122.3 

1.5  /irn 

36 

± 

3 

/irn 

Infragold® 

89.9 

9.4  /irn 

38 

± 

3 

/irn 

Sandblasted  Aluminum 

67.7 

2.3  /irn 

38 

± 

3 

/irn 

Graphite 

37.3 

3.5  /irn 

41 

± 

3 

/irn 

White  paint 

36.8 

1.7  /irn 

41 

± 

3 

/irn 

Spectralon® 

13.8 

Unprofiled 

45 

± 

3 

/irn 

Transmissive  Sample 

White  paint 

56.4 

1.7  /irn 

63  ±  3 

/irn 

Iteration  (count) 

Figure  3.  Enhancement  comparisons  of  specular  and  diffuse  regions  of  reflection.  En¬ 
hancement  is  plotted  after  each  optimization  with  the  CCD  placed  in  the  specular 
region  at  45°  from  the  reflector  surface  normal  (blue).  Specular  is  defined  as  a  reflec¬ 
tion  angle  of  45°  ±  0.5°  from  the  surface  normal.  Diffuse  is  defined  as  15°  off  specular. 
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The  brushed  aluminium  reflector  produced  the  highest  level  of  enhancement  the 
initial  background  speckle  intensity.  All  of  the  diffuse  reflectors  used  in  the  origi¬ 
nal  experiment  showed  at  least  an  order  of  magnitude  enhancement  in  the  focused 
spot  for  N  =  1,  024.  The  original  experiments  also  demonstrated,  using  the  iter¬ 
ative  optimization  algorithm,  that  the  focused  spot  in  the  observation  plane  could 
be  positioned  anywhere  within  the  CCD,  displaced  up  to  8.3  mm  diagonally,  at  a 
distance  of  50  ±  0.5  cm  from  the  reflector,  and  experienced  less  than  a  13.5%  drop 
in  peak  intensity  of  the  focused  spot.  When  the  diagonal  distance  is  increased  to 
108.2  mm,  the  peak  intensity  of  the  spot  is  decreased  by  58.6%.  Using  these  diagonal 
displacements,  the  specular  region  was  defined  with  an  angle  of  incidence  and  angle 
of  reflection  equal  to  45°  ±  0.5°  and  the  diffuse  region  was  defined  as  15°  off  specular. 
Figure  3  shows  the  enhancement  achieved  by  each  material  for  the  specular  and  dif¬ 
fuse  regions.  The  final  optimized  SLM  for  Figure  3  contained  N  =  256  segments,  but 
with  pre-optimization  the  total  process  performed  340  segment  optimizations.  Figure 
3  also  shows  that  highly  reflective  materials  that  are  strongly  scattering  (i.e.  less 
specular)  maintained  enhancement  better  at  larger  displacements  of  the  focus  spot 
than  did  the  more  specular  samples  [9]. 

The  speckle  pattern  produced  by  a  fixed  SLM  phase  map  and  diffusely  reflecting 
sample  varied  with  time.  The  decorrelation  of  the  speckle  pattern  is  a  function  of  both 
the  material  properties  of  the  diffuse  reflector  and  the  laboratory  test  conditions.  The 
correlation  coefficient  was  plotted  with  time  for  the  brushed  aluminum,  Spectralon®, 
and  transmissive  white  paint  samples,  and  is  shown  in  Figure  4.  As  expected,  the 
metal  reflectors  maintained  a  much  higher  correlation  coefficient.  Thus,  most  of  the 
decorrelation  for  the  metal  reflectors  is  attributed  to  the  changing  environmental 
conditions  of  the  laboratory,  such  as  temperature,  air  currents  from  the  laboratory 
ventilation  system,  minute  mechanical  vibrations,  and  device  measurement  error  and 
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noise. 


Figure  4.  The  correlation  coefficient  of  the  speckle  pattern  produced  by  an  unopti- 
mized  phase  map  of  as  a  function  of  time.  The  unoptimized  phase  map  applied  to  the 
SLM  consisted  of  N  equal-sized  segments  of  random  phase  values  arranged  in  a  square. 
The  correlation  coefficient  for  N  =  256  and  N  =  1024  are  shown  for  the  brushed  alu¬ 
minum,  transmissive  white-paint-on-glass,  and  Spectralon®  samples.  Reprinted  with 
permission  [9]. 


Previous  work  in  the  transmissive  inverse  diffusion  showed  enhancements  (■ rj )  rang¬ 
ing  from  50  to  1,000  using  both  iterative  and  transmission-matrix  methods  [6,14,19, 
21,34],  The  cause  of  this  wide  range  of  enhancement  values  is  still  currently  being 
investigated,  with  some  of  the  disparity  likely  attributed  to  noise  [34],  The  volumet¬ 
ric  scatterers,  such  as  Spectralon®  and  white-paint  samples,  showed  a  much  higher 
degree  of  decorrelation  compared  to  the  metals.  This  decorrelation  over  time  adds  a 
degree  of  uncertainty  to  the  process,  where  previously  optimized  SLM  segments  are 
in  fact  suboptimal  by  the  time  the  process  is  complete.  These  suboptimal  segments 
cause  a  decrease  in  the  overall  spot  enhancement  [9]. 
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2.4  Summary 


Path-length  difference  is  the  main  cause  of  the  light  scatter  and  resultant  speckle 
intensity  patterns  in  both  transmissive  and  reflective  inverse  diffusion.  In  the  trans¬ 
missive  case  the  path-length  difference  is  caused  by  the  large  number  of  random 
paths  through  the  turbid  medium.  In  the  reflective  case,  the  path  length  differences 
are  caused  by  the  wavelength-sized  surface  height  fluctuations  of  the  material.  Ul¬ 
timately,  path-length  differences  are  expressed  as  phase  differences  in  the  light  field 
and  can  therefore  be  compensated  for  through  phase  modulation. 

The  proof-of-concept  experiments  for  reflective  inverse  diffusion  provided  no  math¬ 
ematical  models  or  simulations.  The  experiment  assumed  the  ideal  enhancement 
performance  was  comparable  if  not  equal  to  the  transmission  case.  A  diffraction- 
based  model  is  developed  in  Chapter  III  and  surface  statistics  are  used  to  develop  an 
expression  for  the  ideal  enhancement  for  reflective  inverse  diffusion. 
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III.  Mathematical  Model  for  Reflective  Inverse  Diffusion 


In  this  chapter,  rough  surface  models  with  Gaussian-distributed  surface  height 
fluctuations  are  used  to  derive  an  expression  for  ideal  enhancement  for  reflective  in¬ 
verse  diffusion.  The  original  reflective  inverse  diffusion  experiments  took  over  10  hours 
to  optimize  the  SLM  phase  map  of  N  =  1,  024  segments  so  that  it  refocused  light  to 
a  single  CCD  segment  [9].  A  mathematical  model  to  drive  numerical  simulations  was 
needed  to  further  develop  reflective  inverse  diffusion.  Such  a  model  is  also  developed 
in  this  chapter. 


3.1  Enhancement  for  Reflective  Inverse  Diffusion 


Light  scattering,  whether  by  transmission  through  or  reflection  from  a  medium, 
is  a  linear  process.  Despite  complexity  and  unknown  material  properties,  it  can 
be  considered  deterministic  as  long  as  the  medium  is  static  [33].  Therefore,  the 
relationship  developed  for  transmissive  inverse  diffusion  for  the  observed  total  field  in 
the  target  area  is  also  valid  for  the  reflective  case.  The  field  at  the  mth  position  was 
given  in  Equation  (1)  [5],  repeated  here  as 


N 

Em  'y  ^  bnn  An( 

n=  1 


Jfin 


(1) 


where  tmn  represents  the  mnth  element  of  the  transmission/reflection  matrix,  and 
Ane^n  is  the  field  from  the  nth  segment  of  the  SLM  used  to  perform  the  phase 
modulation.  The  source  field  is  segmented  by  the  SLM  into  N  segments  and  reflected 
off  the  scattering  medium.  The  field  in  the  observation  plane  consists  of  a  linear 
combination  of  these  N  segments  of  the  source  field  at  each  of  the  M  observation 
positions  [5].  Normalizing  Equation  (1)  to  intensity,  An  =  G=,  and  the  observed 
intensity  at  the  rnth  position  was  given  by  Equation  (2),  also  repeated  here, 
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n—  1 


(2) 


Intensity  enhancement  was  defined  by  Vcllekoop  for  transmissive  inverse  diffusion 
in  Equation  (3)  as  [7], 


V  = 


(I opt) 


(3) 


\Jrnd) 

which  is  a  simple  ratio  of  the  average  intensity  of  the  optimized  spot  divided  by  the 
average  intensity  of  the  unoptimized  random  speckle  [6].  For  transmission  through 
a  random  medium,  the  tmn  terms  follow  a  statistically  independent  complex  circular 
Gaussian  distribution  [13]  with  properties  that  can  be  used  to  simplify  Equation  (3) 
and  express  the  ideal  intensity  enhancement  as  a  function  of  the  total  number  of  SLM 
segments  [7],  given  in  Equation  (4)  as, 


Tltranmissive  ^  (-^  1)  T  1  ~  4^’  (^) 

The  multiple  random  paths  light  travels  when  transmitted  through  a  scatter¬ 
ing  medium  are  not  present  in  the  reflection  model.  This  lack  of  complex  circular 
Gaussian  statistics  produces  a  different  expression  for  ideal  enhancement.  Using  a 
simplified  geometric  approximation,  the  surface  properties  are  modeled  as  a  constant 
average  surface  reflectivity  and  phase  delay  that  is  related  to  the  reflector  surface 
height  fluctuations  [10].  Using  this  model,  Equation  (2)  is  rewritten  for  the  reflection 
case  as 

2 

(5) 

where  r  is  the  average  surface  reflectivity  and  6,mn  is  the  phase  delay  caused  by 
the  surface  height  fluctuations  of  the  material.  Rough  surfaces  have  surface  height 


Im~  N 


N 

r  ej61" 

n=  1 
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fluctuations  spanning  multiple  wavelengths  [10].  A  rough  surface  modeled  with  a 
Gaussian  surface  height  distribution  with  a  standard  deviation  greater  than  half  a 
wavelength,  will  produce  a  uniformly  distributed  phase  when  wrapped  to  the  interval 
[— 7r,  7r].  Thus,  the  phase  term  of  the  reflector  material,  9mn,  is  assumed  to  have  a 
uniform  distribution  over  [— n,  ti\.  The  intensity  achieves  a  maximum  value  when  the 
SLM  phase  delay  is  set  to  cancel  the  phase  delay  imposed  by  the  surface  height  of 
the  reflector,  (f>n  =  —9mn,  for  a  given  m,  then 


(I max) 


2n, 


(6) 


where  the  angled  brackets,  again,  denote  the  ensemble  average. 

The  unoptimized  random  speckle  background  can  be  expressed  as  a  fixed-length 
random  phasor  sum, 


(And)  —  (  jy- 


N 

E 

n= 1 


r  e 


jdr, 


=  r 


N 


EpjOmn 
n=  i 


(7) 


which  can  be  written  as  a  complex  number  with  amplitude  A  and  phase  0, 


(And)  =  r2  ^|Ae?e|2^  =  r2  (A2)  ,  (8) 

where  A  and  0  are  both  random  variables  and  (A2)  is  the  second  moment.  The 
probability  density  function  of  A  is  given  by  [10], 

Pa(<x)  =  4vr2A  pJq  A(2vrap)dp,  (9) 

where  a  >  0  is  the  magnitude  of  the  random  phasor  sum,  p  =  from  the  2- 
dimensional  joint  characteristic  function  of  the  random  phasor  sum  [10],  and  J$  is 
the  zero-order  Bessel  function  of  the  first  kind.  Using  Equation  (9),  Mathematica® 
approximated  the  second  moment,  (A2)  =  1.  Substituting  the  results  from  Equa- 
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tions  (6)  and  (8)  into  Equation  (3)  shows  the  ideal  enhancement  for  reflective  inverse 
diffusion  is  equal  to  the  total  number  of  SLM  segments  N, 

(Jmax)  AT  /  ln\ 

Vreflective  /  j  \  (10) 

\J-rnd / 

3.2  Diffraction  Theory 

3.2.1  Rayleigh-Sommerfeld  Diffraction. 

The  first  Rayleigh-Sommerfeld  diffraction  formula  can  be  seen  as  a  mathematical 
representation  of  the  Huygens-Fresnel  principle.  The  formula  shows  that  the  held  at 
an  aperture  of  extent  E,  U (Pi),  can  be  described  as  an  infinite  number  of  point  sources 
each  producing  its  own  spherical  wave.  The  observed  held,  U(Pq),  is  the  sum  of  the 
contributions  from  these  point  sources.  The  integral  form  of  the  Rayleigh-Sommerfeld 
diffraction  formula  is  given  by, 

u(P0)  =  ^JJu(Plf^^comis,  HD 

s 

where  A  is  the  optical  wavelength,  k  is  the  propagation  constant  and  equals  The 
integration  is  over  the  enclosed  surface  S  that  contains  P\ .  The  source  and  observation 
planes  are  separated  by  a  screen  that  is  opaque  outside  of  the  aperture  represented 
by  area  E.  The  distance  between  the  points  in  the  source  and  observation  planes  is 
represented  by  r0 1,  and  6  is  the  angle  between  the  surface  normal  of  the  aperture  and 
the  ray  to  the  observation  point  [35,36]. 

Consider  a  system  of  two  parallel  planes,  the  source  plane  with  coordinates  (£,77)5 
and  the  observation  plane  with  coordinates  (x,y).  Since  6  is  the  angle  between 
the  source  plane’s  surface  normal  (n)  and  the  vector  between  Pq  and  P\  (f*oi),  then 
cos (6)  =  — ,  where  the  2-axis  and  n  are  colinear.  Thus,  the  held  at  any  point  in  the 
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observation  plane  (x,  y )  can  be  thought  of  as  a  sum  of  contributions  from  an  infinite 
number  of  point  sources  in  the  source  plane  where  the  limits  of  integration  would 
depend  on  the  aperture  [35,36]. 


U2(x,y) 


z 
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Uifov) 


exp  (j  Ax 01 ) 


d£  dr/ 


01 


(12) 


rra  =  r-  -I  ic  ■  Cr  +  (y  -  nr  (13) 

Rayleigh- Sommerfeld  diffraction  requires  only  the  distance  between  the  points  in 
the  source  and  observation  planes  be  much  greater  than  the  wavelength  (roi  3>  A). 
This  is  the  most  accurate  scalar  wave  optics  diffraction  formula,  and  is  the  basis  for 
Fresnel  and  Fraunhofer  diffraction  [36]. 

3.2.2  Fresnel  Diffraction. 

The  square  root  in  the  distance  term  roi  in  Equation  (13),  can  cause  problems 
with  analytic  solutions  to  diffraction.  The  square  root  can  be  eliminated  by  using  the 
first  two  terms  of  its  binomial  expansion  to  approximate  r0i. 


VT+b  =  l  +  l-b  -  \b2  + ... 

z  o 


(14) 


roi 


=  ^Z2  +  (x  -  t)2  +  (y  -  T])2  =  z^l+  +  (^—^j 


(15) 


r oi  ~  2 


1 1 1  i 1 

2  l  z  2 


z 


(16) 


The  roi  term  appears  in  both  the  exponent  and  denominator  of  Equation  (12).  Since 
the  denominator  is  a  squared  term,  it  is  sufficient  to  reduce  the  approximation  of 
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r0i  ~  z.  The  r0i  in  the  exponent  is  more  susceptible  to  small  changes,  thus,  Equation 
(16)  is  used  [35],  to  produce 


2  rr  exp(jfc  [l  +  l(!rf)2  +  l(t2)2l) 

U2(x,y)  =j-xjj  lh K.t,) - i - 1 - ’-U-dtidr,.  (17) 


ejkz 

Now  simplify  by  moving  — —  outside  the  integral  to  get 

zz 

U2(x,y)  =  JJ  Ui(£,  rj)  exp  [(a;  -  O'2  +  (y  ~  v)2]^j  dr/.  (18) 

By  expanding  the  quadratic  in  the  exponential  and  moving  outside  the 

integral,  Equation  (18)  can  be  re-written  as, 


^jkz 

U2(X,y)  =  e—e^W) 


jXz 


Ui(Z,y)e^ 


ji-jew) 


d^  dr/.  (19) 


Equation  (19)  is  the  Fresnel  diffraction  integral.  The  second  exponential  in  the 
integral  is  the  Fourier  transform  kernel.  Given  the  approximation  on  roi,  this  shows 
that  the  field  in  the  observation  plane  is  proportional  to  the  two-dimensional  Fourier 
transform  of  field  in  the  source  plane,  £A(£,  r/),  multiplied  by  a  quadratic  phase  factor. 

Determining  the  accuracy  of  the  Fresnel  integral  is  accomplished  by  looking  at 
the  approximation  of  rm.  The  accepted  condition  for  accuracy  is  when  the  b2 / 8  term 
that  was  dropped  from  the  roi  approximation  causes  a  maximum  phase  change  much 
less  than  1  radian.  The  accuracy  condition  is  met  when  the  distance  between  the 
source  and  observation  plane,  z,  satisfies  Equation  (20)  [35], 


z3^^[(x~02  +  (y~v)2} 


2"l  2 
max ' 


(20) 


As  an  example,  for  the  propagation  from  the  diffuse  reflector  to  the  CCD,  the 
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diffraction-limited  spot  produced  by  the  SLM  and  positive  lens  (Li)  represents  the 
source  aperture.  The  radius  of  the  diffraction  limited  spot  produced  by  the  SLM  and 
positive  lens  is  given  by  [35], 


Q  = 


A/ 


(21) 


where  /  is  the  focal  length  of  the  lens,  and  dseg  is  the  width  of  the  SLM  segment. 
For  N  —  1,  024  SLM  segments,  then  dseg  =  240  fin i  and  the  radius  of  the  diffraction 
limited  spot  is  approximately  1.3  mm.  The  dimension  of  the  observation  plane  is  the 
6.9  mm  by  4.3  mm  of  the  CCD.  The  wavelength  of  the  HeNe  laser  source  used  in 
the  experiment  is  approximately  633  nm.  Applying  Equation  (20)  to  the  parameters 
of  the  original  experiment,  the  Fresnel  diffraction  approximation  would  require  a 
distance  z  >>  12  cm.  The  observation  plane  is  said  to  be  in  the  near  field  when  the  z 
satisfies  Equation  (20). 


3.2.3  Fraunhofer  Diffraction. 

Further  assume  that  the  distance  between  the  source  and  observation  plane  is 
large  enough  that  the  quadratic  exponent  of  the  first  exponential  in  Equation  (19) 
can  be  approximated  to  be  zero,  so  that 


z>  +  (22) 

Then,  the  Fresnel  diffraction  integral  in  Equation  (19)  can  be  re-written  as  the  Fraun¬ 
hofer  diffraction  integral, 


U2(x,y)  =  %-*&*+*) 


j\z 


£M£,7?)e-^(a*+OT)  d£  drf. 


(23) 


The  distance  requirement  imposed  by  Equation  (22)  can  be  very  large.  Due  to  the 
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large  distances  required  for  Fraunhofer  diffraction,  it  is  commonly  referred  to  as  the 
far  held  approximation  [36].  Using  the  original  experimental  parameters,  the  radius 
of  the  diffraction  limit  spot  is  1.3  mm  as  the  source  aperture  with  a  HeNe  laser 
wavelength  of  633  nm;  Equation  (22)  requires  a  distance  of  z  >>  17  m  for  Fraunhofer 
diffraction  approximation  to  be  valid. 

It  is  possible  to  reduce  the  distance  requirement  of  Equation  (22)  by  the  use  of  a 
positive  lens.  The  positive  lens  of  focal  length  /,  imparts  a  quadratic  phase,  given  by 
Equation  (24),  onto  the  source  held  directly  in  front  of  the  lens, 


tA{£,v)  =  exP 


-j^(e+v2) 


The  total  source  held  is  then  given  by 


(24) 


U0(tv)  =  U^V)  P((,ri)  (25) 


P(Z,V) 


{1,  inside  the  lens  aperature. 
0,  otherwise. 


(26) 


where  P(£,  rf)  is  the  pupil  function  of  the  lens  aperture.  Substituting  Equation  (25) 
into  the  Fresnel  diffraction  integral  in  the  Equation  (19)  gives, 


U2(x,y) 


ajkz 


jXz 


ej£-z(x2+y2) 


£/i(£,  rj)  P(Z,rj)  tA(£,  d£d rj 

(27) 


where  the  quadratic  phase  imparted  by  the  lens  cancels  out  the  quadratic  phase  in 
the  Fresnel  diffraction  integral  when  z  =  f.  This  leaves  only  the  Fourier  transform  of 
the  source  held  and  pupil  function  product.  If  the  extent  of  the  lens  pupil  function 
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P  is  much  larger  than  that  of  the  source  field  Ui,  then  the  effect  of  the  lens  aperture 
can  be  neglected  leaving  only  the  source  field, 


U2(x,y)  =  J f  Ui(Z,  v)  e~j^+yr,)  d£d rj.  (28) 

Equations  (28)  and  (23)  are  identical  for  z  —  /,  but  Equation  (28)  did  not  require 
use  the  stringent  requirements  on  z  of  Equation  (22)  [35]. 

The  Fresnel  number  provides  a  more  straight  forward  method  to  determine  when  it 
is  appropriate  to  apply  either  the  Fresnel  or  Fraunhofer  approximation  to  diffraction. 
Let  w  represent  the  half  width  of  the  aperture,  the  Fresnel  number  (Np)  is  then  given 
by 


ic2 

7VF=-.  (29) 

For  a  Fresnel  number  approximately  equal  to  one  (Np  ~  1),  it  is  common  to  apply 
Fresnel  Diffraction.  In  some  cases,  Fresnel  diffraction  may  be  reasonably  accurate  for 
Fresnel  numbers  as  high  as  30.  The  Fresnel  number  for  Fraunhofer  diffraction  can  be 
determined  from  Equations  (22)  and  (29), 


2  » 


Ke + v 


2N|  7TU72 

max  /l  ' ' 


1  »  — —  =  nNF. 
\z 


(30) 


The  factor  of  n  is  often  dropped  and  the  Fresnel  number  for  Fraunhofer  diffraction 
is  accepted  as  simply  much  less  than  one  (Np  1).  Outside  these  regions  Raylcigh- 
Sommerfeld  diffraction  should  be  used  to  maintain  accuracy  [36]. 
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3.2.4  Transfer  Functions. 


Diffraction  integrals  can  become  very  complex,  very  quickly.  A  different  approach 
is  to  treat  the  held  propagation  through  the  system  as  a  linear  spatial  filter.  The 
impulse  response  of  Rayleigh- Sommerfeld  diffraction  is  given  by  [35], 


h(x,y)  =  - 


z  exp  (jkr) 


3  A 


(31) 


where  z  is  the  propagation  distance,  A  is  the  wavelength,  the  wavenumber  is  k  =  -y, 
and  r  =  z2  +  x2  +  y2.  The  impulse  response  can  be  used  to  rewrite  the  Rayleigh- 
Sommerfeld  diffraction  integral  in  Equation  (12)  as  a  convolution, 


uout(x,  y)  =  jxJJ  Uin &  h(x  -£>y-  v)  d^dri-  (32) 

Applying  the  Fourier  transform  convolution  theorem,  Equation  (32)  can  be  rewritten 
as, 


Uout(x,y)  =  j^Jr~1{r{Uin(£,ri)}r{h(x,y)}} 

=  ^F-1{Jr{Uin(ti,v)}H(fx,fy)} 


(33) 


where  the  transfer  function,  H(fx,  fy),  is  the  Fourier  transform  of  the  impulse  response 
h(x,y)  [35], 


3.3  Mathematical  Model 

Scalar  diffraction  theory  was  used  to  develop  a  mathematical  description  of  reflec¬ 
tive  inverse  diffusion.  The  propagation  distance  of  50  cm  from  the  diffuse  rehector  to 
the  CCD  is  greater  than  the  near-held  distance  of  12  cm  calculated  in  Section  3.2.2 
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using  Equation  (20).  Whether  this  satisfies  the  much  greater  than  requirement  of 
Equation  (20)  is  open  for  interpretation.  However,  since  the  number  of  operations  is 
identical  in  MATLAB®  regardless  of  which  transfer  function  is  used,  the  more  accu¬ 
rate  Rayleigh- Sommerfeld  diffraction  formula  is  used  for  propagation.  The  transfer 
function  for  Rayleigh- Sommerfeld  diffraction  is  [35], 

H(fx,  fy)  =  exp  (^jkz^/l-  (A/a;)2  -  ( Xfy )2^  ,  (34) 

where  z  is  the  propagation  distance,  A  is  the  optical  wavelength,  k  is  the  wavenumber 
27t/A,  and  fx  and  fy  are  the  respective  horizontal  and  vertical  spatial  frequencies  of 
the  source  field.  Using  the  transfer  function  to  propagate  the  source  field  Usrc(x,y )  a 
distance  z,  the  observed  field  is  given  by  [36] 

Uobs{x,y)  =  j^Jr~1{Jr{Usrc{x,y)}  H{fx,fy)},  (35) 

where  (x,  y )  are  the  Cartesian  coordinates  orthogonal  to  ^-direction  corresponding  to 
the  spatial  frequencies  (fx,  fy). 

The  field  at  each  pixel  of  the  SLM  when  illuminated  by  a  plane  wave  can  be 
considered  simply  as  the  phase  delay  applied  at  that  pixel.  Using  Equation  (35),  the 
field  is  propagated  from  the  SLM  to  just  prior  to  lens  Li  in  Figure  2,  normalized  with 
constant  phase  terms  removed  is 

UE1(x,y)  =  F-1{F{UsLM(x,y)}H(fx,fy)  |*=Zl},  (36) 

where  Z\  is  the  distance  from  the  SLM  to  lens  L\.  The  Fourier  transform  property  of 
lens  L\  is  then  used  to  determine  the  field  at  the  diffusely  reflecting  sample  located 
at  the  lens  focus, 
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U7ample(X,y)  =  eJ ^{U^X,  y)} 

=  ej^u'W)T{  T-\T{UsLM(x,y)}  H{fx,  fy) \z=Zl}}  (37) 

=  ej^u2W) F{USlm(x,  y)}  H(fx,  fy) \z=Zl. 

Lens  L\  causes  a  coordinate  transformation  to  the  (u,  v )  plane  which  is  related  to  the 
spatial  frequencies  by  u  =  —  and  v  —  — ,  where  z  —  Z\  [35]. 

The  held  at  the  sample  given  by  Equation  (36)  is  multiplied  by  e^0^u’v')  which  rep¬ 
resents  the  phase  scattering  properties  of  the  rehector.  The  result  is  then  propogated 
to  the  CCD  detector  in  the  observation  plane  using  Equation  (35), 


UCcd(u,v)  = 


(38) 


where  (x,y)  are  the  coordinate  axes  of  the  SLM  and  (fx,  fy)  are  the  respective  hori¬ 
zontal  and  vertical  spatial  frequencies  of  the  held  at  the  SLM.  The  coordinate  axes  of 
the  rehective  sample  are  (u,v),  and  (fu,  fv)  are  the  respective  horizontal  and  vertical 
spatial  frequencies  of  the  held  at  the  rehector.  The  focal  length  of  lens  L\  is  /,  and  Z\ 
and  Z-2  are  the  distances  from  the  SLM  to  the  lens  and  from  the  rehective  sample  to 
the  CCD,  respectively.  This  transform  relationship  is  unique  to  the  rehective  inverse 
diffusion  setup  in  Figure  2.  Transmissive  inverse  diffusion  uses  microscope  objectives 
on  both  sides  of  the  scattering  sample  that  beam-contracts  the  light  from  the  SLM 
onto  the  sample  and  then  re-images  the  light  leaving  the  sample  onto  the  CCD  (see 
Figure  2  in  reference  [6]).  The  differences  in  the  experimental  setups  create  differ¬ 
ent  transform  relationships  for  rehective  and  transmissive  inverse  diffusion.  However, 
since  both  processes  are  linear,  Equation  (1)  is  valid  for  both. 
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3.4  Sample  Plane  Pixel  Size 


The  sample  spacing  of  the  simulated  SLM  was  fixed  at  15  /mi  to  match  the  pixel 
spacing  of  the  physical  SLM.  The  pixel  spacing  establishes  the  maximum  spatial 
frequency  of  [/5]max  =  ^  =  33-3fWT-  The  result  for  Uvlmax  is  identical  since  the 
pixel  spacing  is  uniform  both  vertically  and  horizontally.  The  spatial  frequency  step 
size  is  A/e  =  -A,  where  L g  is  the  length  of  the  source  plane  in  the  ^-direction.  Since 
the  source  plane  is  square,  =  Lv  =  Lv,  and  A/t  =  A fv.  The  size  Li  must  at 
a  minimum  be  equal  to  the  dimension  of  the  SLM,  7.68  mm;  however  it  is  often 
recommended  that  L i  be  two  or  three  times  this  size  [36].  LIsing  the  values  of  [f^]max 
and  A/e,  the  range  of  spatial  frequencies  is  given  by, 


1  .  1  .  1  1 

^  [  2A£  ’  Li  2A£  Li 

with  a  similar  result  for  fv.  Thus,  the  spatial  frequency  resolution  can  be  increased 

by  increasing  L^.  This  is  accomplished  by  padding  the  simulated  SLM  input  with 

zeros. 

The  spacing  at  the  diffuse  reflector  is  then  Ax  =  yy-  It  is  possible  to  pad  L\ 
large  enough  to  achieve  wavelength-scale  values  for  Ax  in  an  attempt  to  simulate  the 
sampling  the  imperfections  of  the  reflector  surface.  This,  however,  should  be  avoided 
as  scalar  diffraction  theory  would  no  longer  be  a  valid  approximation,  and  the  increase 
in  the  size  of  L\  causes  a  substantial  increase  in  processing  time. 
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3.5  Simulations 

A  diffraction-based  propagation  model  allowed  for  examination  of  the  field  at 
any  point  in  the  simulation.  Using  MATLAB®  to  calculate  the  field  at  the  CCD 
using  Equation  (38),  an  Intel  Core  i7®  computer  with  8  gigabytes  of  RAM  could 


complete  over  2  optimizations  per  second,  which  was  more  that  64  times  faster  than 
experimental  methods  with  available  equipment.  The  MATLAB®  source  code  for  the 
simulations  can  be  found  in  Appendix  A.  This  model  was  validated  using  a  mirror 
as  a  test  case. 

The  diffusely  reflecting  sample  was  replaced  by  a  mirror  to  validate  the  propagation- 
based  simulation.  In  simulation,  the  mirror  is  considered  a  perfect  reflector,  which 
eliminates  the  ej0(-x,y>  term  from  Equation  (38).  Qualitatively,  with  a  mirror  posi¬ 
tioned  at  the  focus  of  the  positive  lens,  the  creation  of  a  focused  spot  on  the  CCD 
simply  requires  shifting  the  focus  of  the  positive  lens  the  distance  from  the  mirror  to 
the  CCD.  This  is  accomplished  by  applying  a  negative  lens  phase  screen  to  the  SLM. 
The  phase  screen  for  a  lens  was  given  in  Equation  (24).  Using  geometric  optics,  the 
focus  of  the  negative  lens  is  given  by, 

f SLM  =  1 - 1 - Zi,  (40) 

7  -  7+z~2 

where  Z\  is  the  distance  from  the  SLM  to  the  positive  lens,  /  is  the  focal  length  of 
the  positive  lens,  and  Z2  is  the  distance  from  the  mirror  to  the  CCD. 

The  spot  produced  by  the  mirror  test  case  was  captured  by  the  CCD  and  compared 
with  the  spot  simulated  using  Equation  (38).  Figure  5  shows  both  the  measured  and 
simulated  spots  and  includes  the  center  horizontal  cross  sections  of  each.  There  was  a 
20- jum  difference  between  the  FWHM  diameter  of  the  measured  and  simulated  spots. 
The  difference  is  attributed  to  uncertainties  in  the  experimental  distances  between 
the  SLM,  lens,  reflector,  and  CCD  that  result  a  small  defocus  error  in  the  measured 
spot. 

Computer  simulations  of  the  original  reflective  inverse  diffusion  experiments  cor¬ 
roborate  the  experimental  data.  The  experimental  enhancement  is  plotted  per  it¬ 
eration  of  the  algorithm  in  Figure  6(a).  The  effects  of  speckle  decorrelation  were 
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Figure  5.  Simulation  Validation.  The  diffusely  reflecting  sample  was  replaced  by  a 
mirror  to  compare  measured  and  simulated  intensity  patterns.  A  negative  lens  phase 
screen  was  used  to  position  the  focus  onto  the  CCD. 
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Figure  6.  Enhancement  performance  comparison  of  simulated  and  experimental  data. 
The  N=256  iterations  are  highlighted  in  gray.  The  N=l,024  iterations  are  highlighted 
in  yellow,  (a)  Enhancement  achieved  in  laboratory  experiments  for  each  of  the  six 
reflector  materials.  Reprinted  with  permission  [9].  (b)  Enhancement  achieved  in  simu¬ 
lations  of  reflective  inverse  diffusion  with  zero-mean  Gaussian  distributed  random  phase 
fluctuations  added  to  the  reflector  model.  With  the  exception  of  a  =  8°,  all  reflector 
models  are  unimodular  with  a  uniformly  distributed  random  phase.  For  a  =  8°,  both 
the  uniformly  distributed  random  phase  model  (red)  and  the  circular  Gaussian  model 
(blue)  are  shown. 
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simulated  by  adding  a  zero-mean  random  Gaussian-distributed  phase  to  the  reflec¬ 
tor  model  each  time  the  intensity  was  calculated.  The  standard  deviation  (a)  of 
the  random  phase  fluctuation  was  dependent  on  the  severity  of  decorrelation  time  of 
the  material.  Brushed  aluminum  maintained  the  highest  correlation  coefficient  over 
time,  thus  experienced  the  lowest  amount  of  enhancement  loss  attributed  to  speckle 
decorrelation.  Simulations  of  the  metal  samples  corresponded  to  a  Gaussian  phase 
fluctuation  with  standard  deviation  of  0.5°  to  2°.  The  more  diffuse  samples,  such 
as  graphite  and  white  paint,  matched  with  a  standard  deviation  of  3°  to  5°.  Figure 
6(b)  shows  enhancement  per  iteration  for  a  given  standard  deviation  of  the  Gaussian 
phase  fluctuation. 

Spectralon®  is  a  bulk  scatter,  and  due  to  a  high  depth  of  light  penetration, 
was  considered  similar  to  the  random  paths  experienced  by  light  when  transmitted 
through  a  scattering  medium.  Simulations  of  both  the  circular  Gaussian  model  from 
transmissive  inverse  diffusion  and  uniform  phase  model  produced  similar  enhance¬ 
ment  with  the  same  phase  fluctuation  standard  deviation  deviation  of  a  =  8°  (see 
Figure  6(b)).  The  enhancement  for  the  uniform  phase  model  does  initially  rise  faster 
than  the  circular  Gaussian  model,  which  is  expected  based  on  the  slopes  of  Equations 
(4)  and  (10).  However,  due  to  the  decorrelation  effects  and  limited  number  of  phase 
values  used  in  the  optimization  process,  neither  model  was  significantly  better  than 
the  other  at  simulating  the  Spectralon®  results. 

3.6  Observation 

The  iterative  algorithm  of  reflective  inverse  diffusion  is  a  simple  brute  force  method 
for  determining  phase  screens  that  refocus  light  after  reflection,  but  it  also  has  several 
drawbacks.  Only  a  small  number  of  phase  values  can  been  tested,  and  each  additional 
phase  value  that  is  added  requires  an  additional  N  intensity  measurements.  Each  time 
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the  target  segment  in  the  observation  plane  changes  locations,  the  entire  algorithm 
must  be  completely  restarted.  The  reflection  matrix  method  discussed  in  the  Chapter 
IV  can  solve  both  of  these  issues  by  measuring  the  phase  contributions  of  the  reflector. 
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IV.  Reflection  Matrices 


4.1  Introduction 

Dual  photography  is  a  method  for  mathematically  interchanging  the  positions  of 
a  camera  and  light  source.  This  is  accomplished  by  taking  advantage  of  Helmholtz 
reciprocity.  A  pixelated  light  source  is  used  to  illuminate  a  scene  while  a  camera, 
with  a  different  viewpoint  than  the  light  source,  records  images  of  light  reflection 
and  scattering  from  the  target  scene.  These  images  are  used  to  construct  a  trans¬ 
port  matrix  that  maps  light  source  pixels  to  camera  pixels.  The  transport  matrix 
is  used  to  construct  an  image  of  the  scene  from  the  perspective  of  the  light  source. 
The  reconstructed  image  contains  information  that  was  not  directly  visible  from  the 
perspective  of  the  camera  [2], 

Target  scene  information  hidden  from  the  camera  but  recovered  using  dual  pho¬ 
tography  must  be  visible  from  the  perspective  of  the  light  source.  This  requirement 
limits  the  usefulness  outside  the  author’s  original  work  of  relighting  scenes.  Indi¬ 
rect  photography  is  a  method,  developed  at  the  AFIT,  that  sought  to  co-locate  the 
light  source  with  the  camera,  by  eliminating  the  line-of-sight  requirement  of  the  light 
source,  while  maintaining  the  ability  to  reconstruct  images  not  directly  visible  from 
the  camera.  The  initial  experiments  using  multiple  reflection  radiometric  models 
achieved  limited  success  [3,4], 

The  basis  for  imaging  with  light  scattered  by  transmission  through,  or  reflect  off, 
an  object  was  established  by  Freund.  The  random  scatterer  was  treated  as  a  complex 
field  that  interfered  with  incident  light  to  produce  a  distinct  speckle  pattern.  It  was 
theorized  that  by  properly  modifying  the  incident  light,  the  scatterer  could  be  used 
to  simulate  various  optical  instruments  such  as  a  lens  or  a  mirror  [37].  In  the  case  of 
the  mirror,  this  would  allow  the  imaging  of  objects  using  diffusely  reflected  light  from 
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a  wall  and  to  effectively  see  around  corners.  Proper  development  of  these  methods 
could  have  profound  impact  in  remote  sensing  and  other  fields. 

Phase  modulation  techniques  have  demonstrated  the  ability  to  shape  a  wavefront, 
causing  light  to  refocus  to  a  single  point  after  transmission  through  a  scattering 
medium.  This  process,  called  “inverse  diffusion”  by  the  authors,  was  also  capable 
of  controlling  the  location  of  the  focused  spot  in  the  observation  plane  by  adjusting 
the  wavefront  shape.  The  scattering  of  light  is  a  linear  process,  whether  caused  by 
transmission  through  a  medium  or  reflection  off  its  surface  [6,7].  Applying  the  concept 
of  inverse  diffusion  to  indirect  photography,  refocusing  light  from  the  first  diffuse 
reflection  has  the  potential  to  simplify  the  method  back  to  that  of  dual  photography. 
From  the  observation  plane  perspective,  the  origin  of  the  light  source  is  the  first 
diffuse  reflector,  thus  satisfying  the  line  of  sight  requirement  of  dual  photography.  It 
has  been  demonstrated  that  the  iterative  techniques  for  wavefront  shaping  developed 
by  Vellekoop  et  al.  can  be  adapted  to  work  in  reflection  [38]. 

Matrix  methods  for  inverse  diffusion  have  been  developed  with  significant  ad¬ 
vantages  over  iterative  techniques  for  refocusing  light  after  transmission  through  a 
turbid  medium  [14,15,17,19,21,23].  The  number  of  required  intensity  measurements 
are  significantly  reduced  compared  to  iterative  methods  and  the  phase  information 
collected  allows  refocusing  of  the  light  to  any  location  in  the  observation  plane,  in¬ 
cluding  producing  multiple  foci  simultaneously.  Adapting  these  methods  for  reflection 
would  potentially  bridge  the  gap  between  dual  photography  and  indirect  photogra¬ 
phy  allowing  images  to  be  produced  from  reflectively  scattered  light  and  see  around 
corners. 
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4.2  Methodology 


4.2.1  Parallel  Wavefront  Optimization. 

The  method  used  to  measure  the  RM  was  based  on  work  by  Yoon  et  al.  for 
measuring  transmission  matrices  of  turbid  media  [17].  Yoon’s  method  was  based 
on  parallel  wavefront  optimization  method  by  Cui  and  expanded  to  measure  the 
entire  transmission  matrix,  rather  than  optimize  to  a  single  point  [17,19].  Parallel 
wavefront  optimization  uses  interference  between  reference  and  signal  fields  produced 
by  the  SLM  to  extract  the  phase  information  of  the  RM  [17,19].  The  reference  field  is 
generated  by  static  segments  of  the  SLM  and  the  signal  field  generated  by  modulated 
segments  of  the  SLM. 

The  SLM  is  divided  into  N  segments  and  the  segments  are  separated  into  two 
interleaved  groups.  Each  segment  of  Group  1  is  modulated  at  a  distinct  frequency, 
while  all  the  segments  in  Group  2  are  held  static  as  shown  in  Figure  7(a)  (similar 
to  Figure  2  from  reference  [17],  modified  for  clarity).  However,  assuming  the  SLM 
is  illuminated  with  a  plane  wave  and  all  the  segments  are  initially  set  to  zero  phase 
delay,  using  Equation  (1),  the  field  at  the  mth  position  becomes, 

N  G 

E  =  f  pLV  (AD 

J—/m  /  J  °mn  i  /  J  ump'~'  i  / 

n=G+ 1  p=  1 

where  G  is  the  total  number  of  segments  in  Group  1,  and  (j)p  =  copt  is  the  phase  value 
of  the  SLM  segments  at  time  t.  The  first  summation  in  Equation  (41)  represents  the 
static  reference  field  produced  by  the  segments  in  Group  2  and  simplifies  to  a  complex 
number  that  will  be  represented  by  Rna ■  The  intensity  at  the  rnth  position  is  then 
given  by, 
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(a) 


Figure  7.  SLM  segment  diagram  for  parallel  wavefront  optimization,  (a)  Group  1 
segments  (black)  are  modulated,  each  at  a  distinct  frequency,  while  Group  2  segements 
(white)  are  static,  (b)  Group  2  segments  (white)  are  modulated,  each  at  a  distinct 
frequency,  while  Group  1  segments  (black)  are  static,  (c)  All  Group  1  segments  (black) 
are  modulated  at  u>i,  while  all  Group  2  segments  (white)  are  modulated  at  0J2  (similar 
to  Figure  2  from  reference  [17],  modified  for  clarity). 
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The  desired  matrix  information  is  contained  in  the  second  term  of  Equation  (42). 
Since  the  desired  information  occurs  at  specific  frequencies,  the  R*m2tmp  coefficients 
can  be  determined  using  the  Fourier  transform  of  the  intensity  measurements  [17]. 
The  frequencies  for  the  modulated  SLM  segments  are  given  by, 


cop 


G  +  p 
AG 


^ S  ) 


(43) 


where  cos  is  the  sampling  frequency.  This  ensures  that  up  G  (\ujs,  |u;s]  and  the 
harmonic  frequencies  (c oq  —  cop)  G  [—  \us,  \ojs]  do  not  overlap  [17].  This  requires 
AG  intensity  measurements  for  proper  resolution  in  the  frequency  domain  in  order  to 
extract  all  the  R*m2tmp  terms.  The  roles  of  the  SLM  segment  groups  are  then  reversed, 
Group  1  segments  are  held  static  while  Group  2  segments  are  modulated  as  shown 
in  Figure  7(b).  An  additional  4 (IV  —  G)  intensity  measurements  are  recorded  making 
the  total  number  of  intensity  measurements  41V  [17].  Fourier  transforms  are  used 
to  extract  the  R*mltmp  terms  from  the  second  set  of  intensity  measurements.  This 
results  in  an  M  x  N  matrix  of  R*mxtmp  values  that  contains  the  phase  information 
relating  each  segment  of  the  SLM  to  each  segment  of  the  CCD  detector.  The  matrix 
coefficients  are  simply  amplitude  scaling  and  phase  shifts  to  provide  the  proper  linear 
combination  of  the  fields  from  the  SLM  segments  that  produces  the  field  at  the  CCD 
in  the  observation  plane.  Thus,  the  coefficients  are  by  definition  unitless. 
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4.2.2  Reference  Field  Phase  Matching. 


The  matrix  of  values  collected  consists  of  R*n2tmp  terms  from  the  optimization  of 
Group  1  segments  and  R*nltmp  terms  from  the  optimization  of  Group  2  segments. 
There  is  also  an  ambient  static  reference  field,  Rm 0,  produced  by  real  world  device 
and  laboratory  imperfections.  Thus,  the  matrix  values  can  be  represented  by  [17], 


Group  1  .  (RmO  T  Rrnz)  (bnl)tm2,  *  '  •  Rmc) 

Group  2  :  (Rm o  +  Rml)*  (tm(G+l):  Ln(G+2),  •  •  •  )  tmN  )  • 


(44) 


The  different  reference  fields  provide  a  different  amplitude  scaling  and  phase  shift  to 
the  two  SLM  segment  groups.  Since  the  SLM  is  a  phase-only  device,  the  different 
amplitude  scaling  is  of  little  concern.  However,  the  phase  difference  between  the  two 
reference  fields  must  be  determined  in  order  to  build  the  final  RM. 

The  phase  shift  between  the  two  reference  fields  can  be  determined  by  modulating 
both  SLM  segment  groups  at  different  frequencies.  This  is  similar  to  the  previous  step, 
except  all  the  SLM  segments  of  a  given  group  are  modulated  at  the  same  frequency 
as  shown  in  Figure  7(c).  The  two  frequencies  used  for  this  part  of  the  process  can  be 
determined  using  Equation  (43)  with  G  =  2.  The  field  produced  at  the  CCD  is  then 
a  sum  of  the  ambient  field,  Rm 0,  and  the  modulated  reference  fields  R,rn\  and  Rm 2. 
The  equation  for  intensity  becomes  [17], 


Im  —  \Em\2  —  |  RmO  +  Rml +  Rm2e^2t\ 

2 

=  \Rmk\2  +  RrnoR^e-^  +  Rm0R*m2e-i“2t  +  R^R*^-^-^  (45) 
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Fourier  transforms  are  used  to  extract  the  reference  field  coefficients  from  Equation 
(45).  Eight  intensity  measurements  are  required  for  proper  frequency  domain  reso¬ 
lution.  These  coefficients  are  used  to  determine  the  phase  shift  necessary  to  bring 
the  optimized  Group  2  SLM  segments  in  phase  with  the  optimized  Group  1  SLM 
segments.  The  phase  shift  for  the  optimized  Group  2  SLM  segments  is  given  by  [17], 


G  —  {Rm 0  T  Rm 2)  ( RmO  T  Rml) 


—  pWn  0  ~r 
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Using  the  coefficients  extracted  from  Equation  (45),  the  phase  shift  term  ^  is  calcu¬ 
lated  and  applied  to  the  optimized  Group  2  SLM  segments. 


Group  1  : 
Group  2 : 
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The  term  is  now  common  to  both  groups.  This  phase  only  term  shows 

that  the  reference  held  phase  contribution  to  the  RM  is  simply  a  global  phase  shift 
and  can  be  ignored  [17].  The  RM  values  are 


Group  1  .  |7?m0  T  77m2|  (tmh  ^m2>  >  tmG) 

Group  2  .  |77m0  T  7?m^|  (7m(G+ip  7m(G+2p  j  tmNji 


(48) 


which  shows  the  reference  fields  for  the  two  optimized  groups  still  have  different 
amplitude  scalings,  but  since  the  SLM  is  phase  modulation  only,  this  difference  is 
ignored.  The  argument  of  Equation  (48)  provides  the  phase  information  for  the 
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RM  [17]. 


SLM  Channel  Index 


Figure  8.  The  RM  has  dimensions  M  x  N  and  contains  the  phase  information  of  the 
light  at  each  CCD  channel  from  each  SLM  channel.  To  refocus  the  light  to  a  specific 
CCD  channel  the  corresponding  row  of  the  RM  values  is  conjugated  to  bring  all  the 
incident  light  on  that  CCD  channel  in  to  phase.  The  conjugate  RM  row  vector  is 
reshaped  according  to  the  dimensions  of  the  SLM  and  applied. 


The  RM  is  constructed  from  Equation  (48)  as  show  in  Figure  8.  This  matrix 
describes  the  effect  the  scattering  sample  has  on  the  phase  of  the  incident  light  from 
every  SLM  channel  to  every  observation  plane  channel.  To  focus  light  into  a  specific 
channel,  extract  the  row  corresponding  to  the  desired  observation  plane  channel  and 
reshape  the  row  vector  to  a  matrix  corresponding  to  the  SLM  pixels  (as  shown  in 
Figure  8).  These  phase  values  are  then  conjugated  and  applied  to  the  SLM  [17]. 

The  RM  is  capable  of  generating  multiple  foci  simultaneously.  This  is  done  by  a 
linear  combination  of  phases  from  multiple  rows  of  the  reflection  matrix.  The  phase 
values  applied  to  the  SLM  are  given  by, 

K  i 

SLMfoa  =  -=  Utmpy  ,  (49) 

m=  1  v 

where  f  indicates  the  complex  conjugate,  and  Ztmp  is  the  argument  of  the  complex 
value  tmp.  Thus,  the  phase  applied  to  each  SLM  segment  is  an  average  of  the  values 
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from  the  RM  for  the  individual  observation  plane  channels.  Different  foci  can  be 
given  different  weights  by  replacing  the  -^=  with  am,  where  =  1-  The  SLM 

phase  map  for  multiple  foci  with  different  scaling  factors  is  then  [17], 

I< 

SLM foci  =  ^2  (Ltmp)]  .  (50) 

m=l 

Since  magnitude  information  of  the  RM  is  not  available,  and  only  phase  modulation 
is  used  to  refocus  light,  different  observation  plane  channels  will  experience  different 
levels  of  enhancement.  Equation  (50)  provides  a  method  of  controlling  the  relative 
enhancement  of  each  of  the  targeted  observation  plane  channels  [17]. 

4.2.3  Simulation  Setup. 

Simulations  were  conducted  in  MATLAB®  using  Equation  (38)  to  determine  the 
held  and  intensity  in  the  observation  plane.  The  MATLAB®  source  code  is  available 
in  Appendix  A.  The  distances  Z\  and  Z2  from  Figure  2  are  15  cm  and  50  cm  respec¬ 
tively  and  the  lens  focus  is  /  =  500  mm.  A  single  SLM  pixel  is  15  /mi  square  and  is 
represented  by  a  single  pixel  in  the  simulation.  The  512  by  512  array  of  SLM  values 
is  then  zero-padded  to  provide  the  “guard-area”  for  diffraction  to  prevent  artifacts  at 
the  edges  of  the  simulated  SLM  [36].  Due  to  the  spatial  Fourier  transform  produced 
by  the  positive  lens,  the  amount  of  zero-padding  directly  affects  the  spatial  frequency 
resolution  which  determines  the  pixel  size  in  the  observation  plane.  The  pixel  size  in 
the  observation  plane  is  given  by, 


A  xobs  =  —  (51) 

where  L  is  the  length  of  the  zero-padded  SLM.  For  all  simulations,  256  zeros  where 
padded  to  all  sides  of  the  simulated  SLM  for  a  total  of  1024  x  1024  pixels,  each  15- /mr 
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square  for  a  total  area  of  15.36  mm  square.  This  established  the  observation  plane 
pixel  size  at  20.6  /irn.  The  pixel  size  in  the  observation  plane  is  of  particular  signifi¬ 
cance  because  the  observation-plane  and  sample-plane  spacings  are  the  same.  Thus, 
the  pixel  size  in  the  observation  plane  represents  the  minimum  lateral  correlation 
length  of  the  simulated  scattering  sample. 

4.2.4  Experimental  Setup. 

The  experimental  setup  was  shown  in  Figure  2.  A  Thorlabs  5-mW  HeNe  laser 
(A  =  632.8  nm)  with  a  linearly  polarized  output  is  used  as  the  illumination  source. 
The  Meadowlark  P512  is  a  reflective  SLM  that  is  7.68  mm  by  7.68  mm,  consisting  of 
512  by  512  pixels,  each  capable  of  over  16,000  discrete  phase  levels  over  a  2n  phase 
stroke.  Feedback  is  provided  by  a  Thorlabs  4070-GE  CCD  with  2048  by  2048  pixels 
and  a  pixel  pitch  of  7.5  /irn  by  7.5  /irn. 

The  laser  is  first  expanded  and  collimated  to  fill  the  SLM.  A  non-polarizing  beam 
splitter  (NPBS)  is  used  to  maintain  normal  incidence  with  the  SLM.  The  distance 
Z\ ,  from  the  SLM  to  lens  L\  is  15  ±0.5  cm.  After  modulation,  the  beam  is  focused  by 
a  500-mm  lens,  L\1  onto  the  scattering  sample  placed  at  the  lens  focus.  The  sample 
is  a  6-inch  square  of  polished  aluminum.  The  CCD  was  placed  50  ±  0.5  cm  from  the 
scattering  sample.  Both  SLM  and  CCD  are  controlled  through  MATLAB®.  The 
RM  values  are  determined  using  MATLAB®  fast  Fourier  transform  (FFT)  of  the 
intensity  measurements  from  the  CCD. 

4.3  Results  and  Analysis 

4.3.1  Simulation  Results. 

The  simulated  speckle  patterns  did  not  extend  across  the  entire  1024  by  1024 
pixel  area  of  the  observation  plane.  Since  these  simulations  involve  large  arrays  of 
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Figure  9.  Observation  plane  intensity  patterns  (a)-(c)  are  computer  simulations  and 
(d)-(e)  are  laboratory  results  of  the  focal-plane  system  in  Figure  2  with  the  polished 
almuninum  sample.  Above  each  spot  are  the  (row,  column)  coordinates  of  the  given 
observation  plane  segment.  Phase  maps  generated  by  the  RM  are  used  to  refocus  light 
to  single  or  multiple  segments  in  the  observation  plane,  (a)  A  simulated  single-segment 
enhancement  of  77  =  348  over  background  speckle,  (b)  Simulation  of  two  foci  optimized 
simultaneously  at  (72,  72)  and  (96,  96)  with  enhancements  of  77  =  47  and  56,  respectively, 
(c)  Simulation  of  three  foci  generated  with  increased  background  speckle  at  (72,72), 
(96,96),  and  (120,120)  with  enhancement  values  of  77  =  26,  32,  and  29,  respectively,  (d) 
Demonstrated  single-segment  enhancement  of  77  =  18  over  background  speckle,  (e)  Two 
segments  are  optimized  simultaneously,  measured  enhancement  for  both  foci  is  77  =  10. 
(f)  Three  foci  are  generated  with  increased  background  speckle  at  (24,24),  (32,32),  and 
(40,40),  enhancement  values  are  77  =  7,  6,  and  7,  respectively. 
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data,  only  the  center  768  by  768  pixels  of  the  observation  plane,  where  the  speckle 
patterns  are  the  strongest,  were  used  for  determining  the  RM.  To  further  reduce  the 
computer  memory  requirements  of  the  simulations,  the  windowed  observation  plane 
measurements  were  binned  into  M  channels.  Ideal  enhancement  is  given  in  Equation 
(10)  and  depends  on  N,  the  total  number  of  SLM  segments.  The  simulations  show 
a  relationship  to  M ;  in  general,  larger  values  for  M  increase  the  enhancement  of  the 
spot  the  RM  is  able  to  produce.  Since  the  total  number  of  intensity  measurements 
is  (41V  +  8),  any  increase  in  N  increases  the  runtime  by  a  factor  of  4.  However, 
increasing  M  only  increases  the  amount  of  computer  memory  required. 

It  was  observed  that  if  M  was  of  insufficient  size,  then  the  reference  field  phase 
matching  procedure  described  in  section  4.2.2  failed  to  increase  enhancement  of  the 
spot  produced  by  the  RM.  Smaller  M  values  produce  larger  individual  channels.  This 
allows  for  Group  1  and  Group  2  SLM  pixels  to  optimize  to  different  locations  within 
the  larger  observation  plane  channel  and  produce  a  larger  blurred-out  spot  with  lower 
overall  enhancement.  Segmenting  the  768  by  768  pixels  of  the  observation  plane  into 
M  =  36,  864  channels  provided  the  necessary  resolution  such  that  the  reference  field 
phase  matching  process  significantly  improved  spot  enhancement  produced  by  the 
RM.  The  M  =  36,  864  channels  are  arranged  192  by  192  over  the  observation  plane. 
Each  channel  consists  of  4  by  4  pixels  that  represents  a  simulated  area  of  82.4  fin i 
square. 

The  RM  was  used  to  refocus  the  incident  beam  to  1,024  of  the  36,864  observation 
plane  channels.  The  selected  channels  were  evenly  spaced  throughout  the  observation 
plane.  The  average  enhancement  of  the  sampled  channels  was  rjavg  =  307  with  a 
maximum  enhancement  of  rfmax  =  398.  This  represents  roughly  30%-40%  of  the 
maximum  ideal  enhancement  of  Equation  (10).  Simulation  time  was  approximately 
300  seconds  on  an  Intel®  i7  desktop  computer  with  8  gigabytes  of  memory.  Each 
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value  of  m  has  a  (row,  column)  coordinate  corresponding  to  its  location  in  the  192  by 
192  channel  array  representing  the  observation  plane.  The  simulated  RM  was  used 
to  optimize  the  intensity  at  center  channel  (96,96),  m=18,336.  This  single-channel 
optimization,  shown  in  Figure  9(a),  had  an  enhancement  of  rj  =  348  over  background 
speckle.  Multiple  foci  were  also  generated  using  Equation  (50)  and  can  be  seen  in 
Figures  9(b)  and  9(c).  As  the  energy  is  split  between  multiple  foci  the  enhancement 
of  each  spot  is  reduced.  For  the  two  foci  in  Figure  9(b),  the  enhancement  for  the  spots 
at  (72,  72)  and  (96,  96)  are  r/  =  47  and  56,  respectively.  The  enhancement  for  the 
three  foci  at  (72,72),  (96,96),  and  (120,120),  shown  in  Figure  9(c),  are  r/  =  26,  32, 
and  29,  respectively.  As  the  enhancement  decreases,  the  background  speckle  became 
more  visible. 

Simulations  of  scattering  samples  with  larger  correlations  lengths  were  examined. 
By  applying  the  same  random  phase  value  to  groups  of  pixels  in  the  simulated  sam¬ 
ple,  lateral  correlation  lengths  that  are  multiples  of  the  minimum  20.6  /irn  can  be 
simulated.  In  general,  the  larger  correlation  lengths  did  improve  enhancement  for 
some  of  the  observation  plane  channels,  but  as  the  correlation  length  increased,  the 
sample  became  more  specular  and  the  observation  plane  area  where  enhancement 
could  be  achieved  became  smaller.  Doubling  the  lateral  correlation  length  to  41.2  /jrn 
produced  the  best  results  with  N  =  1,024  SLM  segments,  and  an  observation  plane 
of  768  by  768  pixels  divided  into  M  =  16,  384  channels.  The  channels  near  the  center 
of  the  observation  plane  achieved  enhancement  of  rj  «  1,000,  which  is  very  near  the 
predicted  ideal  maximum  of  Equation  (10).  The  enhancement  drops  quickly  as  the 
channels  move  away  from  the  center  producing  an  average  enhancement  of  rjavg  =  145 
for  the  1,024  channels  that  were  measured.  Increasing  the  lateral  correlation  length 
beyond  41.2  //rri  did  not  significantly  improve  simulation  results.  In  some  cases  the 
speckle  patterns  became  less  developed  and  unintended  multiple  foci  were  produced 
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that  reduced  overall  enhancement  as  the  simulated  scattering  sample  became  more 
specular. 

4.3.2  Experimental  Results. 

Intensity  measurements  from  the  CCD  are  used  to  determine  the  RM.  These 
measurements  are  sensitive  to  vibrations  and  changing  environmental  conditions  in 
the  laboratory.  Thus,  the  faster  the  acquisition  of  the  intensity  measurements,  the 
better  the  RM  performs.  The  feedback  from  the  CCD  was  limited  to  the  center  512 
by  512  pixels,  which  limits  the  observation  area  to  3.84  mm  by  3.84  mm.  Limiting  the 
feedback  area  of  the  CCD  both  increased  the  frame  rate  and  decreased  the  memory 
and  computation  time  requirements  to  calculate  the  RM.  The  SLM  was  divided  into 
N  =  1024  segments.  The  feedback  area  from  the  CCD  was  segmented  in  M  =  4096 
channels,  each  channel  was  8-by-8  pixels  with  a  super-pixel  area  of  60  /mi  by  60  /mi. 
The  resulting  RM  was  a  matrix  of  4096  by  1024  elements,  which  required  a  total  of 
4,104  intensity  measurements,  with  an  average  run  time  under  300  seconds. 

Single-channel  optimization  performance  was  measured  by  using  the  RM  to  se¬ 
quentially  refocus  the  beam  to  1,024  of  the  available  4,096  CCD  channels.  The  average 
single-channel  enhancement  was  rjavg  =  9  with  a  maximum  enhancement  of  rjmax  =  19 
over  background  speckle.  Both  single-channel  and  simultaneous  multiple  channel  en¬ 
hancement  are  shown  in  Figure  9(d),  9(e),  and  9(f).  Each  spot  location  is  labeled  with 
its  corresponding  (row,  column)  coordinates.  Single-channel  enhancement  at  (32,  32) 
is  shown  in  Figure  9(d)  with  an  enhancement  of  rj  =  18.  Single-channel  enhancement 
performance  is  less  than  2%  of  the  maximum  ideal  enhancement  predicted  by  Equa¬ 
tion  (10)  and  less  than  5%  of  the  performance  in  the  simulations.  Multiple  foci  were 
generated  along  the  diagonal  using  Equation  (50)  and  are  shown  in  Figures  9(e)  and 
9(f).  Both  foci  in  Figure  9(e)  each  have  an  enhancement  of  rj  =  10,  while  in  Figure 
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9(f)  the  center  spot  enhancement  is  6  and  the  two  outer  foci  have  an  enhancement  of 
7  over  background  speckle. 

The  lower-than-anticipated  enhancement  can  be  partially  attributed  to  non-ideal 
attributes  of  the  laboratory  equipment  and  environment  such  as  SLM  and  CCD  fill 
factors,  vibrations,  and  temperature  changes.  The  primary  cause  for  the  degraded 
enhancement  performance  is  attributed  to  the  lateral  correlation  length  of  the  alu¬ 
minum  plate,  which  is  assumed  to  be  much  less  than  the  minimum  simulated  correla¬ 
tion  length  of  20.6  /jrn.  A  sample  with  a  smaller  lateral  correlation  length  produces  a 
higher  degree  of  scattering  that  the  SLM  segments  are  not  able  to  precisely  match  to 
produce  the  expected  enhancement.  Simply  increasing  the  number  of  SLM  segments 
is  not  likely  to  increase  performance.  The  current  experimental  setup  shown  in  Fig¬ 
ure  2.  It  is  the  2-dimensional  spatial  Fourier  transform  of  the  held  at  lens  Li  that 
is  being  scattered  by  the  sample.  The  held  at  the  lens  is  determined  by  propagating 
the  held  from  the  SLM  to  the  lens.  By  increasing  the  number  of  SLM  segments  and 
decreasing  their  size,  the  spatial  frequencies  of  the  held  at  the  lens  are  increased. 
Thus,  the  2-dimensional  spatial  Fourier  transform  at  the  sample  covers  a  larger  area, 
further  increasing  the  scatter  and  potentially  decreasing  overall  enhancement. 

The  spatial  Fourier  transform  of  the  held  at  the  lens  provided  a  linear,  predictable, 
and  convenient  relationship  for  simulation  purposes,  but  it  does  not  minimize  the  in¬ 
teraction  with  the  scattering  sample.  To  improve  enhancement  using  this  relationship 
would  require  a  large  number  of  SLM  segments  (IV),  but  the  physical  size  of  the  seg¬ 
ments  would  have  to  increase  in  order  to  minimize  the  size  of  the  spatial  Fourier 
transform  at  the  sample.  In  Chapter  V,  rehective  inverse  diffusion  experiments  will 
attempt  to  increase  enhancement  by  minimizing  the  interaction  area  of  the  light  with 
the  scattering  sample. 
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4.4  Conclusion 


Yoon’s  method  for  measuring  transmission  matrices  has  been  demonstrated  to 
work  in  reflection.  The  RM  maintained  the  ability  to  refocus  light  to  any  of  the 
individual  channels  of  the  CCD  or  multiple  channels  simultaneously.  This  method 
for  measuring  the  RM  required  less  than  20%  of  the  intensity  measurements  of  the 
previous  iterative  method  of  Section  2.3.2  that  was  only  capable  of  optimizing  to  a 
single  channel.  The  enhancement  capabilities  of  the  RM  showed  a  dependence  on 
the  physical  dimensions  and  total  number  of  observation-plane  channels  (M)  that 
was  not  mentioned  in  the  case  of  transmission.  The  overall  enhancement  from  the 
RM  in  the  laboratory  experiments  was  significantly  less  than  anticipated  and  is  at¬ 
tributed  to  a  lateral  correlation  length  of  the  scattering  sample  that  is  smaller  than 
that  which  is  able  to  be  phase  matched  by  the  SLM.  Spot  enhancements  of  up 
to  20  times  background  speckle  were  achieved.  Transmissive  inverse  diffusion  ex¬ 
periments  typically  use  optics  to  produce  a  de-magnihed  image  of  the  SLM  held  at 
the  transmissive  sample  and  then  re-magnify  the  scattered  beam  for  imaging  onto  a 
CCD  [5-7,14,15,17,19,21,33].  Chapter  V  will  examine  adapting  this  approach  for  use 
in  reflection.  Applying  a  de-magnihed  image  of  the  SLM  onto  the  scattering  rehector 
could  minimize  the  interaction  area  and  improve  the  enhancement  performance  of 
the  RM.  However,  implementing  this  approach  in  the  application  of  imaging-around- 
corners  could  prove  to  be  less  practical  than  the  focal-plane  approach.  None  the  less, 
improvement  of  the  enhancement  capabilities  of  the  RM  would  provide  a  valuable 
tool  for  imaging  using  rehected  light  and  effectively  seeing  around  corners. 
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V.  Reflection  Matrix  Measurement 


5.1  Introduction 

The  majority  of  surfaces  can  be  considered  rough  when  compared  to  the  wave¬ 
length  of  visible  light.  The  microscopic  surface  height  variations  of  the  rough  surface 
cause  incident  light  to  diffusely  scatter.  The  scattering  surface  can  be  treated  as  a 
complex  field  that  interferes  with  incident  light  to  produce  reflected  speckle  patterns. 
If  the  complex  field  of  the  rough  surface  was  known,  a  properly  modified  incident 
beam  could  be  used  to  eliminate  the  scattering  effects  of  the  rough  surface  and  even 
cause  the  light  to  refocus  after  reflection  [37]. 

Transmissive  inverse  diffusion  used  phase  modulation  to  shape  the  wavefront  of 
incident  light  causing  it  to  refocus  after  transmission  through  a  turbid  media  [5- 
7,19,30].  These  methods  were  also  demonstrated  to  work  in  controlling  reflective 
scatter  [9].  Reflective  inverse  diffusion  is  an  iterative  process  with  intensity  feedback 
from  a  CCD  to  search  for  the  SLM  phase  map  that  produced  the  brightest  target  spot 
in  the  observation  plane.  This  required  tens  of  thousands  of  intensity  measurements 
and  limited  the  phase  map  to  a  small  subset  of  available  values. 

Transmission  matrices  have  been  developed  to  map  the  effect  of  the  complex  field 
of  the  scatterer  on  the  incident  light.  These  matrices  provide  the  phase  information 
required  to  control  the  resultant  light  scatter  [14,15,21,23].  Transmission  matrices 
were  measured  with  microscopic  objectives  and  thin  films  of  turbid  media,  resulting 
in  propagation  distance  of  less  than  1  mm  and  a  observation  plane  field  of  view  of  a 
few  hundred  microns.  The  complex  field  representation  is  not  limited  to  transmissive 
scattering,  but  the  reflective  case  represents  a  103  increase  in  propagation  distance 
and  observation  plane  size.  Initial  work  in  reflective  inverse  diffusion  always  placed 
the  rough  surface  reflector  at  the  focus  of  a  positive  lens  [9] .  This  paper  will  continue 


50 


to  examine  the  enhancement  capabilities  of  the  RM  with  the  rough  surface  reflec¬ 
tor  at  lens  focus.  Additionally,  the  RM  will  be  measured  with  the  lens  producing 
a  demagnihed  image  of  the  the  SLM  at  the  rough  surface  reflector,  similar  to  the 
transmissive  case. 

Previous  work  in  the  transmissive  scattering  case  showed  enhancements  (rj)  rang¬ 
ing  from  50  to  1,000  using  both  iterative  and  tranmission  matrix  methods  [6, 14, 19, 
21,34],  The  cause  of  this  wide  range  of  enhancement  values  is  still  currently  being 
investigated,  with  some  of  the  disparity  likely  attributed  to  noise  [34],  Enhancement 
capabilities  of  the  RM  in  laboratory  experiments  and  diffraction-based  simulations 
are  compared  with  surface  roughness,  correlation  length,  and  slope  to  provide  initial 
insight  into  surface  characteristics  that  affect  enhancement.  The  simulations  will  also 
examine  the  effect  on  enhancement  of  simplified  device  measurement  error  for  the 
SLM  and  CCD  along  minute  mechanical  vibrations  that  cause  microscopic  shifts  in 
the  optical  setup. 

Imaging  around  corners  using  reflectively  scattered  light  has  tremendous  applica¬ 
tion  in  remote  sensing.  Previous  work  in  this  area  has  always  required  the  occluded 
scene  to  be  illuminated  by  a  light  source  either  in  the  scene  or  with  direct  line  of  sight 
of  the  scene  [2,39].  This  presents  an  application  problem  since  access  to  the  occluded 
scene  may  not  be  possible.  The  goal  this  RM  research  is  to  provide  a  method  of 
illuminating  the  occluded  target  scene  without  access  or  direct  line  of  sight. 

5.2  Background 

Transmissive  and  reflective  scattering  are  both  linear,  deterministic  processes  as 
long  as  the  scattering  medium  is  static.  In  both  cases,  the  resultant  scattered  held 
can  be  considered  a  linear  combination  of  the  inputs.  Phase  modulation  using  an 
SLM  divides  up  the  incident  light  into  several  several  individual  segments  each  with 
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a  unique  phase.  Every  sampled  location  in  the  observed  speckle  field  is  a  linear 
combination  of  the  field  from  the  individual  SLM  segments.  The  field  at  the  rrith 
segment  of  the  observered  field  is  given  by  Equation  (1)  and  reprinted  here  [5], 
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where  tmn  is  the  mnth  complex- valued  element  of  the  transmission/reflection  matrix 
relating  the  light  from  the  nth  SLM  segment  to  the  mth  segment  in  the  observation 
plane,  and  Ane1^'1  represents  the  amplitude  and  phase  of  the  light  from  the  nth  SLM 
segment.  Normalizing  Equation  (1)  to  intensity  by  setting  An  =  1/Vn ,  the  observed 
intensity  was  then  given  by  Equation  (2),  also  reprinted  here, 
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Intensity  enhancement  (■ rj )  is  a  measure  of  performance  for  controlling  scattered 
light.  Defined  in  transmissive  experiments  as  the  ratio  of  the  average  intensity  of 
the  optimized  segment,  ( Iopt ),  to  the  average  intensity  of  the  unoptimized  random 
segments,  (Irnd)  [7].  Assuming  the  surface  height  imperfections  of  the  reflector  are 
on  the  order  of  a  wavelength  and  follow  a  Gaussian  distribution,  the  reflector  can  be 
modeled  as  an  average  reflectivity  with  a  uniform  phase  [10].  The  reflector  statistics 
can  be  used  to  show  the  expected  maximum  ideal  enhancement,  given  in  Equation 
(10)  and  reprinted  here,  is  equal  to  N,  the  total  number  of  SLM  segments  [9]. 
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The  method  used  to  measure  the  RM  was  based  on  work  by  Yoon  et  al.  for 
measuring  transmission  matrices  of  turbid  media  [17].  Yoon’s  method  was  based  on 
parallel  wavefront  optimization  method  by  Cui  and  expanded  to  measure  the  entire 
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transmission  matrix,  rather  than  optimize  to  a  single  point  [17,19].  Parallel  wavefront 
optimization  uses  interference  between  reference  and  signal  fields  produced  by  the 
SLM  to  extract  the  phase  information  of  the  RM  [17,19].  The  SLM  segments  are 
divided  into  two  groups,  the  Group  1  segments  are  modulated  at  a  separate  frequency 
to  produce  the  signal  held,  and  the  Group  2  segments  are  static,  producing  the 
reference  held.  The  intensity  of  a  segment  in  the  observation  plane  is  a  combination 
of  the  reference  and  signal  helds  and  was  given  by  Equation  (42),  also  reprinted 
here  [17], 
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where  i?m2  is  the  sum  of  the  static  Group  2  segments  that  produce  the  reference 
held.  The  second  term  of  Equation  (42)  shows  the  matrix  coefficients,  tmp,  occur  at 
discrete  frequencies.  The  (. R*m2tmP )  coefficients  can  be  extracted  using  the  temporal 
Fourier  transform  of  the  segment  intensity.  This  process  produces  half  of  the  matrix 
coefficients,  the  roles  of  the  SLM  segments  are  then  switched  to  capture  the  other 
half  of  the  RM  coefficients.  A  third  optimization  is  performed  to  bring  both  halves 
of  the  RM  into  phase  [17]. 

5.3  Methodology 

5.3.1  Laboratory  Experiments. 

The  primary  equipment  for  the  laboratory  experiment  is  a  laser  source,  an  SLM, 
and  a  CCD  for  feedback.  The  laser  source  is  a  Thorlabs  5-mW  HeNe  laser  (A  =  632.8 
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nm)  with  a  linearly  polarized  output.  The  laser  is  beam  expanded  to  fully  illuminate 
the  SLM.  The  Meadowlark  P512  is  a  LCoS  SLM  consisting  of  512  x  512  pixels,  each 
capable  of  over  16,000  discrete  phase  levels  over  a  total  2n  phase  stroke.  The  SLM 
has  an  effective  area  of  7.68  mm  x  7.68  mm  with  pixel  pitch  of  15  /mi.  Intensity 
feedback  is  recorded  using  a  Thorlabs  4070-GE  monochrome  CCD  detector  array 
with  a  resolution  of  2048  x  2048  pixels  and  a  pixel  pitch  of  7.5/jrn. 

A  total  of  five  rough-surface  samples  were  made  from  1-inch  aluminum  squares. 
Each  sample  received  a  different  surface  preparation,  the  first  was  bead  blasted,  while 
the  remaining  four  were  polished  with  100  grit,  220  grit,  320  grit,  and  600  grit  sand 
paper.  Surface  profiles  of  each  sample  was  measured  using  a  Tencor  stylus  profiler  in 
both  the  x  and  y  directions.  The  surface  roughness,  a,  was  calculated  as  the  standard 
deviation  in  surface  height  in  each  direction  and  then  averaged  together  to  produce  a 
single  roughness  value  per  sample.  The  sample  roughness  and  autocorrelation  length 
is  in  Table  2. 

The  RM  of  each  sample  was  measured  in  three  separate  locations  using  the  optical 
setup  shown  in  Figure  2.  This  is  the  same  optical  setup  used  in  the  original  reflective 
inverse  diffusion  experiments  with  the  rough  surface  reflector  was  placed  at  the  focus 
of  lens  L\  [9].  The  reflection  matrix  measurements  were  repeated  using  the  optical 
setup  shown  in  Figure  10.  The  single  lens  system  produced  a  |-scale  image  of  the 
SLM  on  the  rough-surface  reflector.  The  reflection  matrix  was  then  used  to  refocus 
light  individually  to  each  of  the  M  segments  of  the  observation  plane.  The  average 
and  maximum  intensity  enhancements  for  each  sample  are  also  shown  in  Table  2. 

5.3.2  Simulations. 

The  propagation  model  assumes  the  SLM  is  illuminated  by  a  unit  amplitude  ideal 
plane  wave.  The  individual  SLM  segments  are  represented  by  unimodular  complex 
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Figure  10.  Image  plane  optical  setup  for  reflective  inverse  diffusion.  A  vertically 
polarized  HeNe  laser  is  expanded,  collimated,  and  normally  incident  on  the  SLM. 
The  lens  (L i)  produces  a  demagnifled  image  of  the  phase  modulated  beam  onto  the 
rough  surface  reflector  and  the  reflected  speckle  pattern  is  recorded  by  the  CCD.  The 
mirror  (Mi)  and  the  NPBS  are  used  to  allow  normal  incidence  with  the  SLM.  The 

demagnification  of  the  system  is  —  sa  -. 
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coefficients.  The  phase  of  these  coefficients  represents  the  phase  delay  imparted  by 
the  SLM  segments  onto  the  incident  plane  wave.  Absorption  and  transmission  are 
neglected  when  modeling  the  aluminum  rough-surface  reflectors;  thus,  the  samples 
are  also  modeled  as  unimodular  complex  numbers.  The  random  phase  imparted  by 
the  rough  surface  is  related  to  the  surface  height  fluctuations  or  surface  roughness. 
Assuming  a  Gaussian  surface  height  distribution  with  a  standard  deviation  equal 
or  greater  than  A/2>  ®ie  phase  distribution  becomes  uniform  when  wrapped  to  the 
interval  [— it,  7r]  [10]. 

Diffraction-based  models  were  simulated  using  the  MATLAB®  2-dimensional  fast 
Fourier  transform  (FFT)  and  the  Rayleigh- Sommerfeld  transfer  function  [36].  For 
the  focal-plane  optical  setup  in  Figure  2,  the  complex  field  at  the  observation  plane 
is  calculated  using  the  Rayleigh-Sommerfeld  transfer  function  and  spatial  Fourier 
transform  property  of  the  lens.  The  held  in  the  observation  plane  was  given  by 
Equation  38  and  reprinted  here  [9], 


UCcd(u,v )  = 


(38) 


where  (x,  y )  are  the  coordinates  of  the  source  plane  at  the  SLM  and  (u,  v )  are  the 
coordinates  of  the  sample  plane  at  the  lens  focus.  The  horizontal  and  vertical  spatial 
frequencies  of  the  source  and  sample  planes  are  ( fx,fy )  and  fv),  respectively. 
The  focal  length  of  lens  L\  is  /,  and  z\  and  Z2  are  the  distances  from  the  SLM 
to  the  lens  and  from  the  reflective  sample  to  the  CCD,  respectively.  The  rough 
surface  is  represented  by  e®®® ,  where  the  phase  function  9(u,  v )  is  a  random  uniform 
distribution  from  [ — 7T,  7r] . 

Ideal  imaging  is  assumed  for  the  configuration  in  Figure  10.  Propagation  from  the 
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SLM  to  the  rough  sample  is  reduced  to  a  global  phase  shift  and  can  be  ignored.  The 
rough  surface  then  randomizes  the  phase  of  the  ideal  image,  which  is  then  propagated 
to  the  CCD  in  the  observation  plane  using  the  Rayleigh-Sommerfeld  transfer  function. 
Thus  the  held  in  the  observation  plane  is  given  by, 


UccdM  =  T-l{7{e^USLM(-u,  -v)}H(fuJv) \Z=Z2}  (52) 

Si  Si 

where  s0  is  the  object  distance  from  the  SLM  to  the  lens  L i,  and  s*  is  the  image 
distance  from  the  lens  Li  to  the  rough  surface  sample.  The  MATLAB®  source  code 
for  the  simulations  is  available  in  Appendix  A. 

5.3.3  Segment  Size. 

The  number  of  SLM  segments,  N,  is  directly  proportional  to  the  maximum  ideal 
enhancement  of  the  refocused  spot  on  the  CCD.  The  parallel  wavefront  modulation 
method  used  to  measure  the  RM  requires  4IV  +  8  intensity  measurements  [17].  For 
all  the  experimental  RM  measurements,  the  SLM  was  divided  into  1,024  equal-sized 
segments,  arranged  32  x  32,  which  required  4,104  intensity  measurements  to  measure 
the  RM.  This  was  the  largest  number  of  segments  that  could  be  equally  sized  and 
utilized  the  entire  SLM  area. 

The  segment  size  of  the  CCD  also  affects  enhancement.  Parallel  wavefront  modu¬ 
lation  measures  the  RM  in  two  separate  halves.  If  the  CCD  segment  size  is  too  large, 
the  two  halves  of  the  RM  may  focus  to  separate  locations  within  the  same  CCD 
segment  resulting  in  a  larger  blurred  spot  with  lower  overall  enhancement.  The  only 
indication  of  a  CCD  segment  that  is  too  large  is  the  lack  of  enhancement  increase  dur¬ 
ing  the  third  step  of  parallel  wavefront  optimization,  phase  matching  the  two  halves 
of  the  RM.  Selecting  a  CCD  segment  size  that  is  smaller  than  necessary  increases 
the  total  number  to  segments  and  increases  the  memory  and  processing  time  of  the 
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RM  measurement.  The  best  CCD  segment  size  was  found  to  be  approximately  l/2 
the  size  of  the  predicted  diffraction-limited  spot  of  the  optical  system. 

The  predicted  spot  size  for  the  optical  system  shown  in  Figure  2  was  calculated  in 
the  original  reflective  inverse  diffusion  proof-of-concept  experiments  [9].  The  radius 
of  the  final  spot  was  given  by, 
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z2  DslM 
f  2  y/N’ 


(53) 


where  z2  is  the  distance  from  the  sample  to  the  CCD,  /  is  the  focal  length  of  lens 
Li,  Dslm  is  the  width  of  the  SLM,  and  N  is  the  total  number  of  SLM  segments. 
Using  Equation  (53),  the  CCD  segment  size  for  the  focal-plane  optical  setup  is  120 
/irn,  which  is  16  x  16  pixels  on  the  Thorlabs  40740M  CCD. 

The  predicted  spot  size  for  the  imaging  system  in  Figure  10  is  calculated  as  the 
diffraction-limited  spot  produced  by  an  aperture  the  size  of  the  demagnified  SLM 
image  that  is  applied  to  the  rough-surface  reflector.  The  radius  of  the  final  spot  is 
then  given  by, 


^  \z2  _  s0  \z2 

^ 2  it  Dslm  Dslm 

where  sQ  is  the  object  distance  from  the  SLM  to  lens  Li,  and  Sj  is  the  image  distance 

from  lens  L\  to  the  rough  surface  sample.  Thus  the  segment  size  for  the  imaging 

system  is  317  /irn,  which  is  42  x  42  pixels  on  the  Thorlabs  4070M  CCD. 


5.4  Results  and  Analysis 
5.4.1  Experiments. 

The  RM  was  measured  three  times  along  the  diagonal  of  each  sample,  top  left 
corner,  center,  and  bottom  right  corner,  for  both  the  focal-plane  and  imaging  optical 
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systems.  The  RM  was  then  used  to  focus  the  incident  beam  to  1,024  different  segments 

of  the  CCD  and  record  the  enhancement  of  each  target  segment  over  the  background 

segments.  The  maximum  enhancement  for  each  sample,  r)max ,  is  the  average  peak 

enhancement  from  each  of  the  three  measured  RMs.  The  average  enhancement  for 

each  sample,  rjavg ,  is  the  mean  enhancement  value  over  all  the  segments  from  all  three 

RMs.  The  rjmax  and  r]avg  for  each  of  the  five  samples  are  recorded  in  Table  2. 

Table  2.  Summary  table  for  the  five  aluminum  samples.  Roughness  (er)  is  the  standard 
deviation  of  the  sample  surface  height.  Correlation  (fc)  is  the  autocorrelation  shift  to 
reduce  the  maximum  value  by  e-1.  Slope  (s)  is  the  RMS  surface  slope.  Correlation 
(fA/2)  is  the  lateral  distance  required  to  change  the  surface  height  by  x/2  at  the  RMS 
slope  (s).  The  samples  average  enhancement  ( r)avg )  and  maximum  enhancement  ( r]max ) 
achieved  with  each  optical  setup  are  included. 


Aluminum 

Roughness 

Correlation 

Slope 

Correlation 

Focal  Plane 

Image  Plane 

Samples 

a 

L 

s 

^A/2 

Vavg 

Vmax 

Vavg 

Vjmax 

100  Grit 

2.13  /mi 

54.25  /mi 

0.143 

2.21  /mi 

1.6 

7.2 

2.2 

10.1 

220  Grit 

1.67  /mi 

35.25  /mi 

0.141 

2.24  /mi 

1.5 

5.8 

3.1 

13.3 

Bead  Blasted 

1.04  /irn 

30.25  /m i 

0.105 

3.00  /mi 

2.2 

8.5 

6.1 

22.5 

320  Grit 

0.63  /irn 

13.75  /mi 

0.107 

2.97  /mi 

1.3 

4.5 

4.3 

16.7 

600  Grit 

0.45  /mi 

16.00  /m i 

0.066 

4.79  /m i 

1.3 

6.0 

5.4 

24.2 

Data  trends  in  Figure  11  show  enhancement  decreases  as  surface  roughness  in¬ 
creases  for  the  imaging  system.  The  enhancement  for  the  focal-plane  system  is  almost 
flat  with  the  average  enhancement  just  slight  larger  than  1;  i.e.  the  target  segment 
intensity  is  not  significantly  higher  than  the  background.  The  difference  between  con¬ 
structive  and  destructive  interference  is  a  surface  height  change  of  thus,  the  rate 
at  which  the  surface  height  changes  is  a  better  predictor  for  enhancement.  The  root 
mean  square  (RMS)  surface  slope  indicates  the  rate  of  the  surface  height  fluctuations 
and  was  determined  from  the  profile  data  from  each  sample.  For  discrete  profile  data, 
the  surface  slope  is  given  by  [11], 


s  = 


1 


K  -  1 


%k  %k—l 
%k  %k—l 


(55) 


where  K  is  the  total  number  of  measurements,  Zk  is  the 


surface  height  from  the 
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Figure  11.  The  surface  roughness  measured  by  the  Tencor  profilometer  is  plotted  with 
enhancement.  Linear  trend  lines  were  added.  For  the  focal  plane  system,  maximum 
enhancement  is  shown  in  red,  while  average  enhancement  is  shown  in  magenta.  For  the 
imaging  system,  maximum  enhancement  is  shown  in  blue,  while  average  enhancement 
is  shown  in  cyan. 
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profilometer,  Xk  is  the  lateral  sample  distance,  and  z'  =  — ^ — — -  is  the  average  slope. 

K  Ax 

The  RMS  slopes  of  the  samples  are  included  in  Table  2.  In  general,  the  samples  with 
larger  RMS  slope  showed  lower  enhancement. 

The  correlation  length,  tc  in  Table  2,  is  defined  as  the  shift  required  to  lower  the 
autocorrelation  of  the  surface  height  profile  by  a  factor  of  e^1  of  the  maximum  value 
[11].  The  apparent  inverse  relationship  between  correlation  length  and  enhancement  is 
counter-intuitive,  but  is  not  conclusive  since  the  samples  with  the  longest  correlation 
length  also  have  the  highest  roughness.  None  of  the  samples  showed  Gaussian  surface- 
height  distributions,  and  with  regards  to  reflective  inverse  diffusion,  the  standard  e_1 
definition  used  for  correlation  length  is  overly  optimistic  due  to  the  small  surface- 
height  difference  between  constructive  and  destructive  interference. 

For  reflective  inverse  diffusion  the  proposed  lateral  dimension,  £\/2  in  Table  2,  is 
defined  as  the  distance  required  to  achieve  a  A/2  change  in  surface  height  given  the 
RMS  slope.  Enhancement  is  plotted  with  £\/2  in  Figure  12.  The  imaging  system 
showed  enhancement  increases  with  i\i2-  The  data  is  inconclusive  in  the  focal-plane 
system  due  to  low  enhancement  values  and  virtually  flat  trends. 

5.4.2  Simulations. 

The  RM  produced  by  the  simulation  of  the  focal-plane  system  in  Figure  2  was 
used  to  refocus  the  incident  beam  to  1,024  of  the  16,384  observation  plane  segments. 
The  rough  surface  reflector  was  made  up  of  individual  segments  with  a  uniform  phase 
distribution.  The  phase  of  each  segment  was  independent  which  leads  to  a  delta 
correlated  sample.  The  segment  size  is  a  measure  of  lateral  correlation,  but  due  to 
the  phase  independence  between  each  segment,  it  was  not  directly  comparable  to  ic 
of  the  tested  aluminum  samples. 

The  phase  of  the  simulated  rough  surface  segments  is  a  uniform  distribution  over 
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Figure  12.  The  lateral  surface  dimension  £w2  is  plotted  with  enhancement.  The  focal- 
plane  system  shows  no  change  in  enhancement  with  ^A/2-  However,  the  average  en¬ 
hancement  of  1  indicates  the  target  segment  has  the  same  intensity  as  the  background. 
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[— 7r,  7r].  If  bounded  over  a  single  wavelength,  the  phase  distribution  gives  a  surface 
height  range  from  [— fvf]-  Since  the  simulated  rough  surface  segments  are  indepen¬ 
dent,  the  surface  height  difference  ( Zk  —  Zk~  1)  in  Equation  (55)  for  surface  slope,  is 
a  random  variable  with  the  same  uniform  distribution  and  an  expected  value  of  A/2. 
Assuming  a  mean  slope  of  z'  =  0,  the  RMS  surface  slope  for  the  simulated  rough 
surface  is  given  by, 


where  Ax  is  the  dimension  of  the  simulated  rough  surface  sample  segment.  Given  the 
of  the  surface  slope  in  Equation  (56),  for  the  simulations  i\/2  =  Ax,  the  simulated 
rough  surface  segment  size. 

The  simulated  SLM  was  zero-padded  as  described  in  section  3.4  to  produce  a 
lateral  surface  dimension  ( i\/2 )  of  10.3  //m.  Decreasing  the  lateral  surface  dimen¬ 
sion  below  10.3  /jm  caused  unacceptable  increases  to  processing  time.  Larger  lateral 
surface  dimensions  were  simulated  by  setting  groups  of  pixels  to  the  same  value.  Ulti¬ 
mately  lateral  surface  dimensions  from  10.3  /jrn  to  660  /jrn  were  simulated.  Maximum 
enhancement  was  achieved  with  £\/ 2  =  41  //rn,  which  produced  an  enhancement  of 
Vmax  =  907  or  88%  of  the  predicted  ideal  maximum  from  Equation  (3).  The  best 
average  enhancement  of  the  1,024  simulated  measurements  was  achieved  with  £\/2  = 
82  /irn,  which  produced  an  enhancement  of  r\avg  =  708.  Increasing  i\/2  beyond  82  /jrn 
caused  both  the  maximum  and  average  enhancements  to  decrease. 

In  the  focal-plane  system,  the  individual  SLM  segments  act  as  individual  aper¬ 
tures.  The  spacing  and  phase  of  the  SLM  segments  create  fringe  patterns  in  the 
diffraction  limited  spot  incident  on  the  rough  surface  reflector.  The  sum  of  all  inter¬ 
ferences  from  the  SLM  segment  pairs  produces  a  wavefront  that  conjugates  the  phase 
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imparted  by  the  rough  surface  reflector  to  produce  a  wave  that  converges  in  the  ob¬ 
servation  plane.  Segment  pairs  that  are  furthest  apart  produce  the  narrowest  fringe 
spacing  of  41  pm.  This  represents  the  smallest  lateral  surface  dimension  (fh/2)  that 
the  focal-plane  system  can  conjugate  and  explains  the  peak  maximum  enhancement 
at  41  pm  for  the  focal-plane  system,  see  Figure  13.  Using  the  RM  to  refocus  light 
ensures  that  the  light  at  the  target  segment  is  all  in-phase,  but  it  does  not  guarantee 
it  is  the  only  segment  where  the  light  is  in-phase.  In  the  focal-plane  system,  as  the 
correlation  length  of  the  simulated  sample  increased,  higher-order  fringes  would  be 
incident  on  the  same  simulated  sample  segment  and  remain  in-phase.  This  produced 
additional  foci  in  the  observation  plane  that  would  increase  background  intensity  and 
decrease  enhancement  as  the  correlation  length  was  increased,  see  Figure  13. 

Simulations  of  the  imaging  system  shown  in  Figure  10  assumed  a  1-mm2  ideal 
image  of  the  SLM  projected  onto  the  rough-surface  reflector.  The  SLM  was  modeled 
with  N=l,024  equal-sized  segments,  arranged  32  x  32  across  the  SLM,  each  with  an 
area  of  31.25  pm  x  31.25  pm.  The  lateral  surface  dimension  of  the  simulated  SLM 
image  is  the  dimension  of  a  single  SLM  segment,  £\/2,slm  =  31.25  pm.  The  rough- 
surface  reflector  was  modeled  with  discrete  segments,  all  with  unit  magnitude  and  a 
uniform  phase  distribution.  The  £\/2  of  the  simulated  rough  surface  was  adjusted  by 
varying  the  number  of  segments  used  in  the  model.  Simulations  were  performed  with 
rough-surface  lateral  dimensions  (fx/2)  as  high  as  500  pm,  and  as  low  as  7.8  pm. 

When  the  lateral  surface  dimension  of  the  simulated  rough  surface  was  identical 
to  the  lateral  surface  dimension  of  SLM  image,  i\/2  =  £\/2,slmi  the  simulated  RM 
achieved  a  maximum  enhancement  of  r)max  =  673  with  an  average  enhancement  of 
rjavg  =  310.  This  represents  30%-65%  of  the  performance  predicted  by  Equation  (3). 
This  performance  drops  to  28%-39%  as  the  correlation  length  of  the  rough  surface 
decreases  to  1/a  of  the  correlation  length  of  SLM  image.  The  best  performance  of 
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Lateral  Surface  Dimension  (^/2)  n m 

Figure  13.  Simulated  correlation  length  fw2  vs  enhancement.  The  red  solid  line  shows 
average  enhancement  of  the  focal-plane  system,  where  the  red  dotted  line  shows  max¬ 
imum  enhancement.  The  solid  blue  line  shows  average  enhancement  of  the  imaging 
system,  with  the  blue  dotted  line  representing  maximum  enhancement. 
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40%-73%  was  achieved  when  l\/2  =  4£A/2 ,slm- 

The  simulations  with  the  longer  lateral  surface  dimension  (i\/2)  produced  greater 
enhancement,  both  the  maximum  achieved  by  a  single  observation  plane  segment, 
and  the  average  enhancement  of  the  simulated  1,024  observation-plane  segments.  The 
enhancement  is  plotted  with  lateral  surface  dimension  (i\/2)  for  both  the  focal-plane 
and  imaging  systems  in  Figure  13.  The  focal-plane  system  achieves  much  higher  peak 
enhancement,  but  only  at  very  specific  correlation  lengths  before  it  rapidly  decreases. 
The  imaging  system  out  performs  the  focal-plane  system  at  shorter  values  of  £\/2  in 
both  maximum  and  average  enhancement.  For  longer  values  of  £\/2  the  enhancement 
of  the  imaging  system  remains  stable  and  does  not  decrease. 

In  the  imaging  case,  the  SLM  is  imaged  and  directly  applied  to  the  rough  surface 
reflector.  The  SLM  segments  conjugate  the  phase  changes  imparted  by  the  reflector 
and  produce  a  converging  wavefront.  As  the  lateral  surface  dimension  {£\/-z)  of  the 
sample  increases,  fewer  adjustments  to  the  phase  map  are  required.  Extended  to  the 
perfect  mirror  case,  the  SLM  phase  map  becomes  the  discretized  phase  function  for  a 
positive  lens  with  the  focal  length  in  the  observation  plane.  The  background  intensity 
does  not  increase  with  correlation  length  as  in  the  focal-plane  system.  As  correlation 
length  increases,  diffraction  from  the  SLM  segments  and  phase  quantization  error  of 
the  SLM  becoming  the  limiting  factors  for  enhancement. 

5.4.3  RM  properties. 

The  RM  measured  with  the  focal-plane  and  imaging  systems  both  maintained  the 
ability  to  refocus  light  to  a  single  CCD  segment  or  multiple  segments  simultaneously 
similar  to  their  transmission  matrix  counterparts.  Multiple  segments  are  enhanced 
simultaneously  using  the  same  process  as  the  transmission  case,  by  using  a  linear 
combination  of  the  rows  of  the  RM  [17].  Simulations  of  the  focal-plane  system  show 
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single-segment  enhancement  in  Figure  9(a),  and  multi-segment  enhancements  shown 
in  Figures  9(b)  and  9(c).  The  imaging  system  with  the  600-grit  aluminum  sample 
demonstrates  single-segment  and  simultaneous  multi-segment  enhancements  and  are 
shown  in  Figures  14(d),  14(e),  and  14(f).  Each  spot  location  is  labeled  with  its 
corresponding  (row,  column)  coordinates.  As  the  enhancement  was  split  over  multiple 
segments,  the  background  intensity  becomes  more  visible. 

5.4.4  Predicted  vs  Measured  Enhancement. 

The  best  enhancement  values  from  the  aluminum  samples  were  10  to  24  times  the 
RMS  background  intensity,  which  is  less  than  10%  of  the  average  enhancement  of  300 
to  400  achieved  in  simulation.  Device  error  does  not  account  for  this  discrepancy. 
For  both  the  imaging  and  focal-plane  simulations,  including  a  random  device  error 
up  to  5%  for  both  the  SLM  and  the  CCD  only  decreased  the  average  enhancement 
by  approximately  10%.  The  RM  is  measured  using  the  temporal  FFT  of  intensity 
measurements  from  the  CCD.  This  makes  the  measurement  process  inherently  noise 
resistant  since  random  device  noise  does  not  occur  at  a  specific  frequency. 

The  aluminum  samples  are  assumed  to  be  stable  for  much  longer  than  the  5  min¬ 
utes  it  takes  to  measure  the  RM.  However  micro- vibrations  in  the  optical  setup  cause 
the  incident  light  to  shift  on  the  rough  surface  reflector.  The  oscillation  of  the  RM 
measurement  area  produces  intensity  measurements  that  belong  to  several  different 
RMs.  For  the  imaging  system  with  a  simulated  rough  surface  lateral  dimension  of 
£\/2  =  7.8  fu n,  incorporating  a  random  vertical  and  horizontal  shift  of  the  sample 
of  23.4  /jrn  reduces  average  enhancement  to  r]avg  =  6  and  maximum  enhancement 
to  r)rnax  =  43.  For  the  focal-plane  system  with  a  simulated  rough  surface  lateral  di¬ 
mension  of  i\j2  =  10.3  /im,  a  simulated  vibration  magnitude  of  20.6  /mi  produces 
maximum  enhancement  of  rjmax  =  12  and  an  average  enhancement  r\avg  =  2.3. 
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Figure  14.  Observation  plane  intensity  patterns  (a)-(c)  are  computer  simulations  of  the 
focal-plane  system  and  (d)-(e)  are  laboratory  results  of  the  imaging  system  in  Figure  10 
with  the  600-grit  aluminum  sample.  Above  each  spot  are  the  (row,  column)  coordinates 
of  the  given  observation  plane  segment.  Phase  maps  generated  by  the  RM  are  used  to 
refocus  light  to  single  or  multiple  segments  in  the  observation  plane,  (a)  A  simulated 
single-segment  enhancement  of  77  =  348  over  background  speckle,  (b)  Simulation  of 
two  foci  optimized  simultaneously  at  (72,72)  and  (96,96)  with  enhancements  of  rj  =  47 
and  56,  respectively,  (c)  Simulation  of  three  foci  generated  with  increased  background 
speckle  at  (72,72),  (96,96),  and  (120,120)  with  enhancement  values  of  77  =  26,  32,  and  29, 
respectively,  (d)  Demonstrated  single-segment  enhancement  of  77  =  12  over  background 
speckle,  (e)  Two  segments  are  optimized  simultaneously,  measured  enhancement  for 
both  foci  is  rj  =  6.  (f)  Three  foci  are  generated  with  increased  background  speckle  at 
(18,6),  (22, 16),  and  (7,22),  enhancement  values  are  rj  =  6,  5,  and  5,  respectively. 


The  lateral  surface  dimension  (^w2)  for  the  aluminum  samples  is  much  smaller 
than  can  be  simulated  without  an  exponential  increase  in  processing  time.  However, 
trends  between  enhancement  and  l\/2  from  the  experiments  and  simulations  were 
extended  and  plotted  for  comparison  in  Figure  15.  The  imaging  system  shows  the  best 
alignment  between  experimental  and  simulation  trends  for  maximum  enhancement. 
The  imaging  system  average  enhancement  trend  lines  for  experimental  and  simulation 
data  have  matching  slope  despite  the  increase  in  separation.  The  experimental  data 
for  the  focal-plane  system  showed  very  low  enhancement  with  predominately  flat 
trends  for  both  average  and  maximum  measured  enhancement.  Similar  enhancement 
levels  are  achieved  in  simulation  for  t\/2  =  10.3  /ini  with  the  addition  of  the  random 
vibration  oscillation  of  20.6  /mi.  Simulations  show  enhancement  increases  for  £\/2  > 
10.3  /mi.  This  could  indicate  the  correlation  length  of  the  aluminum  samples  was  too 
small  for  the  focal-plane  system  to  be  effective  in  the  given  laboratory  conditions. 

5.5  Conclusion 

Yoon’s  method  for  measuring  transmission  matrices  was  successfully  implemented 
to  measure  a  reflection  matrix  (RM).  The  method  only  requires  the  optical  system 
to  be  linear.  Both  the  focal-plane  and  the  imaging  systems  produced  RMs  capable 
of  refocusing  light  onto  the  CCD.  The  RM  is  capable  of  enhancing  a  single  CCD 
segment,  or  multiple  segments  simultaneously. 

Simulations  show  the  focal-plane  system  achieving  higher  levels  of  enhancement 
for  a  very  specific  lateral  surface  dimension  (Tw2).  Outside  of  this  optimum  range, 
the  imaging  system  produces  better  enhancement.  Laboratory  data  shows  the  imag¬ 
ing  system  consistently  produced  higher  levels  of  enhancement  than  the  focal-plane 
system  on  all  the  samples  tested.  This  was  expected  as  simulations  predicted  the 
imaging  system  would  achieve  higher  levels  of  enhancement  at  lower  values  of  i\/2. 
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Figure  15.  Simulations  are  conducted  with  a  random  shift  in  the  sample  prior  to  each 
intensity  measurement.  The  solid  lines  are  experimental  data  trends,  dashed  lines  rep¬ 
resent  simulation  data  trends.  The  simulated  focal-plane  system  is  subject  to  a  random 
20.6  /tm  shift,  maximum  enhancement  is  shown  in  red,  while  average  enhancement  is 
shown  in  magenta.  The  simulated  imaging  system  is  subject  to  a  random  23.4  71m  shift, 
maximum  enhancement  is  shown  in  blue,  while  average  enhancement  is  shown  in  cyan. 
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The  enhancement  for  both  the  focal-plane  system  and  imaging  system  was  signif¬ 
icantly  lower  than  predicted  by  Equation  (10)  for  the  ideal  simulated  results.  The 
primary  cause  of  this  was  minute  mechanical  vibrations  in  the  optical  setup.  Al¬ 
though  enhancement  for  both  methods  were  affected,  incorporating  a  simple  random 
vibration  in  the  simulation  model  showed  that  the  imaging  system  was  capable  of 
greater  enhancement  values  in  this  non-ideal  case.  This  is  corroborated  by  the  higher 
enhancement  values  achieved  by  the  imaging  system  in  the  laboratory. 

Further  research  is  required  to  solidify  the  RM  measurement  method.  Improv¬ 
ing  the  vibration  isolation  of  the  optical  system  will  be  key  to  improving  the  RM 
enhancement  performance.  Mathematical  methods  for  combating  the  effects  of  vi¬ 
bration,  such  as  dynamic  data  driven  applications  systems  (DDDAS)  [40],  should 
also  be  explored  as  a  more  cost  effective  solution.  Ideally,  in  the  imaging  case,  the 
incident  light  is  scattered  by  the  same  part  of  the  rough  surface  regardless  of  the 
which  CCD  segment  is  enhanced.  Since  the  rough  surface  is  assumed  to  be  static, 
the  scattering  is  deterministic  [5].  Chapter  VI  will  investigate  methods  of  exploiting 
the  deterministic  process  to  determine  methods  of  reconstructing  missing  portions  of 
the  RM  from  partial  data,  or  expanding  an  RM  to  include  additional  CCD  segment 
locations  without  requiring  the  RM  to  be  re-measured. 
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VI.  Reflection  Matrix  Reconstruction 


Every  row  of  the  reflection  matrix  (RM)  contains  the  information  necessary  to 
produce  a  converging  wave  after  reflection  off  the  rough-surface  reflector.  Assuming 
the  reflector  properties  are  static  and  the  incident  light  always  interacts  with  the  same 
portion  of  the  reflector  surface,  the  scattering  process  is  deterministic,  which  means 
every  row  of  the  RM  also  contains  surface  height  information  about  the  rough-surface 
reflector.  Since  the  light  interacts  with  the  same  surface  area  for  a  given  RM,  the 
phase  information  that  conjugates  the  phase  changes  imposed  by  the  surface  is  the 
same  for  all  the  rows  of  the  RM.  Thus,  a  large  amount  of  the  information  in  the  RM 
is  redundant,  which  can  be  exploited. 

The  imaging  system  that  was  introduced  in  Chapter  V  produced  greater  levels  of 
maximum  and  average  enhancement  at  smaller  lateral  surface  dimensions  (£w2)  than 
the  original  focal  plane  system  of  Chapter  II.  The  enhancement  of  the  imaging  system 
also  remained  stable  as  G/2  was  increased.  The  imaging  system  also  provided  a  direct 
application  of  an  SLM  phase  map  to  the  rough-surface  reflector,  which  simplified 
pattern  analysis  of  the  RM.  Thus,  the  methods  in  this  chapter  have  been  developed 
for  use  with  the  imaging  system.  However,  these  methods  should  be  adaptable  to  the 
focal  plane  system  provided  the  method  is  modified  to  account  for  the  spatial  Fourier 
transform  relationship  between  the  SLM  and  the  rough-surface  reflector. 

6.1  Geometric  Analysis 

Reflective  inverse  diffusion  is  based  on  the  SLM  ability  to  compensate  for  the 
surface  height  fluctuations  of  the  rough-surface  reflector  and  cause  the  light  to  refocus 
at  a  desired  point  in  the  observation  plane.  The  phase  maps  produced  by  the  RM  can 
be  considered  in  two  parts  that  are  added  together  and  phase  wrapped  over  [— tt,  7r], 


72 


The  first  part  of  the  phase  map  produces  the  converging  phase  front  matching  a 
discretized  phase  function  for  a  positive  lens.  The  second  part  of  the  phase  map 
conjugates  the  phase  shifts  caused  by  the  surface  height  fluctuations  of  the  rough- 
surface  reflector.  The  rough  surface  behaves  as  a  flat  smooth  specular  surface  that 
reflects  the  remaining  converging  phase  front  produced  by  the  phase  map  like  a  mirror. 

Steering  the  focused  spot  to  an  adjacent  observation  plane  segment  requires  a 
625-prad  change  to  the  angle  of  reflection  to  shift  the  focused  spot  312.5  /jrn  in  the 
observation  plane  located  0.5  m  from  the  reflector.  Physically  changing  the  incident 
angle  of  the  system  is  impractical  because  the  incident  beam  would  shift  location  on 
the  reflector  and  invalidate  the  RM.  The  path  length  difference  caused  by  the  small 
change  in  incident  and  reflection  angle  can  be  expressed  as  a  linear  phase  shift  and 
added  to  the  existing  phase  map  of  the  RM.  Using  the  small-angle  approximation, 
the  path  length  is  linear  from  0  to  625  nm  across  the  1-mrn  demagnified  SLM  image 
projected  onto  the  rough-surface  reflector.  Thus,  to  shift  the  target  segment  to  an 
adjacent  segment  in  the  observation  plane  requires  a  linear  phase  shift  from  0  to 
almost  2n,  or  approximately  one  full  wavelength. 

6.2  Reflection  Matrix  Analysis 

The  RM  generated  from  ideal  simulations  was  found  to  contain  linear  phase  shifts 
between  the  rows  of  the  RM  representing  adjacent  observation  plane  segments.  The 
phase  shift  was  comparable  to  the  result  from  the  geometric  analysis:  a  linear,  0-to-27r 
phase  shift,  over  the  width  of  the  simulated  SLM.  The  difference  between  two  phase 
maps  from  the  RM  that  each  refocuses  the  light  to  vertically  adjacent  observation 
plane  segments  is  shown  in  Figure  16(a);  for  two  horizontally  adjacent  observation 
plane  segments,  the  phase  map  difference  is  in  Figure  16(b).  Phase  maps  for  x—  and 
2/ — tilt  calculated  using  the  geometric  approximation  are  shown  in  Figures  16(c)  and 
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16(d).  The  phase  maps  are  similar  to  the  first-order  Zernike  polynomials  for  tilt  [41]. 
However,  since  orthogonality  is  not  a  concern,  a  circular  aperture  is  not  required.  This 
phase  map  can  be  applied  to  shift  the  focus  spot  to  an  adjacent  target  segment  in  the 
observation  plane.  The  linear  phase  difference  was  consistent  throughout  the  entire 
RM,  so  it  is  possible  to  reconstruct  a  new  RM  from  a  single  row  from  an  existing  RM. 


mm  mm 


(c)  (d) 

Figure  16.  (a)  Phase  map  difference  from  two  rows  of  the  RM  corresponding  to  verti¬ 
cally  adjacent  observation  plane  segments,  (b)  Phase  map  difference  from  two  rows  of 
the  RM  corresponding  to  horizontally  adjacent  observation  plane  segments,  (c)  Phase 
map  calculated  from  geometric  analysis  to  shift  target  observation  plane  segment  to 
the  adjacent  segment  below,  (d)  Phase  map  calculated  from  geometric  analysis  to  shift 
target  observation  plane  segment  to  the  adjacent  segment  to  the  right. 


A  single  row  from  an  existing  RM  is  a  phase  map  that  refocuses  light  after  reflec- 
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tion  to  a  specific  segment  in  the  observation  plane.  A  new  RM  can  be  constructed 
by  using  the  tilt  phase  maps  in  Figure  16  to  calculate  phase  maps  for  the  adjacent 
observation-plane  segments.  This  is  repeated  until  the  entire  RM  is  reconstructed. 
In  ideal  simulations,  the  enhancement  performance  of  the  new  RM  is  dependent  on 
the  enhancement  achieved  by  the  original  reference  phase  map  used  in  construction. 
The  original  RM  from  an  ideal  simulation  produced  an  average  enhancement  (rjavg)  of 
310  and  a  maximum  enhancement  (■ r]max )  of  677.  A  reference  phase  map  was  selected 
from  this  RM  that  produces  an  enhancement  (rj)  of  304.  From  this  phase  map,  a 
new  RM  was  constructed  that  produced  rjavg  =  306  and  rjmax  =  546,  approximately 
a  20%  decrease  in  maximum  enhancement.  Using  a  reference  phase  map  that  pro¬ 
duced  an  enhancement  of  658,  the  reconstructed  RM  performance  was  r]avg  =  503  and 

Vmax  —  776,  representing  62%  and  15%  increases,  respectively.  The  RMS  difference 

7 r 

between  the  new  RM  and  the  reconstructed  RM  was  greater  than  —  for  all  the  con¬ 
ducted  ideal  simulations.  This  would  indicate  that  there  exist  many  RMs  for  a  given 
random  rough  surface  reflector  and  the  parallel  wavefront  modulation  method  used 
for  measuring  the  RM  does  not  guarantee  the  best  enhancement  performance. 

6.3  Vibrational  Noise 

Vibrations  in  the  optical  setup  severely  limited  the  enhancement  performance 
of  the  RMs  measured  in  the  laboratory.  These  vibrations  obscured  linear  phase 
tilt  differences  between  phase  maps  for  adjacent  CCD  segments.  Adding  a  random 
vibrational  shift  of  30  //rri  to  the  simulations  made  the  linear  phase  tilt  difference 
that  was  present  in  the  ideal  simulations  all  but  indistinguishable.  However,  even 
under  these  non-ideal  circumstances,  the  method  for  reconstructing  a  new  RM  from 
a  single  phase  map  is  still  valid.  The  original  RM  with  a  30-/im  vibration  produced 
rjavg  —  4  and  rjmax  =  20.  Using  a  phase  map  that  produced  an  enhancement  of  20,  the 
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new  rescontructed  RM  performance  was  increased  to  r]avg  =  7  and  rjmax  =  21,  which 
represents  a  75%  increase  in  average  enhancement. 

6.4  Conclusion 

The  imaging  system  developed  for  measuring  a  RM  directly  applied  the  phase  of 
the  SLM  to  the  rough  surface  reflector  which  allowed  new  phase  maps  to  be  calculated 
by  applying  a  linear  phase  tilt  to  shift  the  target  segment  in  the  observation  plane.  A 
similar  process  should  be  possible  with  the  focal-plane  system  by  applying  the  inverse 
spatial  Fourier  transform  of  the  tilt  phase  maps  in  Figure  16.  However,  this  was  not 
explored  due  to  the  superior  performance  of  the  imaging  system  and  time  constraints. 
Still,  it  should  be  explored  as  future  work. 

Ideal  simulations  show  that  reconstructed  RM  enhancement  performance  is  de¬ 
pendent  on  the  enhancement  of  the  reference  phase  map  used  in  the  construction. 
It  was  demonstrated  that  in  some  cases,  the  reconstructed  RM  can  out-perform  the 
original.  The  RMS  phase  difference  between  the  original  and  reconstructed  RM  indi¬ 
cates  that  multiple  RM  may  exist  for  a  given  surface.  The  increased  performance  of 
the  reconstructed  RM  indicates  that  parallel  wavefront  modulation  method  does  not 
guarantee  maximum  enhancement  performance.  Again,  future  work  could  investigate 
methods  which  guarantee  maximum  enhancement  performance. 

Vibration  in  the  optical  setup  prevents  the  identification  of  the  linear  phase  differ¬ 
ence  between  phase  maps  for  adjacent  CCD  segments.  Introducing  vibration  into  the 
simulations  obscured  the  linear  phase  difference  present  in  the  ideal  simulations.  The 
reconstruction  method  was  still  valid  in  this  case  and  the  reconstructed  RM  showed 
a  substantial  improvment  in  average  enhancement  when  the  reference  phase-map  en¬ 
hancement  was  near  'r]max  for  the  original  RM. 

An  imaging  system  configuration  for  reflective  inverse  diffusion  eliminates  the  need 
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to  measure  the  entire  RM.  However,  measuring  the  RM  with  parallel  wavefront  mod¬ 
ulation  is  currently  faster  than  previous  single-point  iterative  optimization  methods. 
Future  research  should  investigate  more  sophisticated  intensity-based  phase-retrieval 
methods  (e.g.,  steepest  descent,  Gerchberg-Saxton,  conjugate-gradient  method)  for 
single-segment  optimizations.  The  phase  maps  for  the  remaining  segments  in  the 
observation  plane  can  be  determined  using  phase  tilts. 
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VII.  Conclusion 


The  goal  of  this  research  was  to  develop  a  method  of  controlled  illumination  us¬ 
ing  reflected  light  that  could  be  incorporated  into  the  existing  indirect-photography 
configuration  to  effectively  reduce  it  to  the  much  simpler  dual-photography  config¬ 
uration  to  recover  images  around  corners.  The  original  reflective  inverse  diffusion 
research  was  a  simple  proof-of-concept  experiment,  with  little  mathematical  back¬ 
ground,  and  conducted  under  the  assumption  that  the  process  operated  similarly  to 
transmissive  inverse  diffusion  [38].  The  analysis  accomplished  in  this  dissertation 
provided  a  mathematical  basis  for  maximum  expected  enhancement  for  reflective  in¬ 
verse  diffusion  (see  Chapter  III).  Diffraction  theory  led  to  the  development  of  a 
fast-Fourier-transform-based  simulation  that  can  produce  accurate  enhancement  re¬ 
sults  that  take  only  minutes  to  complete  versus  the  original  experimental  runtime  of 
10  hours  (see  Chapter  III).  This  provides  the  means  for  efficient  development  of  new 
reflective  inverse  diffusion  techniques. 

Reflection  matrices  (RMs)  were  developed  in  Chapter  IV  for  reflective  inverse  dif¬ 
fusion  based  on  parallel  wavefront  modulation  [17,19].  Initial  attempts  at  measuring 
a  RM  from  a  rough-surface  reflector  demonstrated  limitations  of  the  original  focal- 
plane  optical  setup  shown  in  Figure  2.  The  process  was  moderately  successful  for  the 
relatively  smooth,  specular,  polished  aluminum  sample.  RMs  provide  the  necessary 
control  of  the  reflected  light  for  use  as  an  illumination  source  in  a  dual  photography 
system. 

The  initial  RM  measurements  led  to  the  development  of  the  imaging-system  op¬ 
tical  setup  in  Figure  10  (see  Chapter  II).  This  setup  is  similar  to  that  used  in  trans¬ 
missive  inverse  diffusion,  on  which  reflective  inverse  diffusion  was  modeled,  but  had 
been  considered  less  practical  in  application  for  the  reflective  case.  The  enhancement 
performance  of  the  two  optical  systems  was  compared  against  surface  properties  of 
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roughness,  correlation  length,  and  surface  slope.  Surface  slope  was  determined  to  be 
the  most  significant  parameter  in  regard  to  enhancement.  The  surface  slope  was  used 
to  define  the  lateral  surface  dimension  (tb/2)  for  reflective  inverse  diffusion.  The  RM 
measurement  time  of  approximately  5  minutes  led  to  the  discovery  of  minute  vibra¬ 
tions  in  the  optical  system  that  significantly  reduced  the  enhancement  performance  in 
laboratory  experiments.  When  small  random  shifts  were  added  to  the  rough-surface 
to  model  these  vibrations,  the  simulation  results  aligned  with  laboratory  experiments. 
The  simulation  capabilities  have  been  continuously  developed  in  both  speed  and  ac¬ 
curacy.  Low  resolution  simulations  of  the  imaging  system  can  now  complete  a  RM 
measurement  in  as  little  as  30  seconds. 

Finally,  in  Chapter  VI,  analysis  of  the  RM  measured  by  the  imaging  system  re¬ 
vealed  a  linear  phase  tilt  relationship  between  phase  maps  that  target  adjacent  seg¬ 
ments  in  the  observation  plane.  The  vibration  noise  that  limited  enhancement  of 
the  laboratory  experiments  also  obscures  these  linear  phase-pattern  differences  of  the 
RMs  measured  in  the  laboratory.  This  relationship  can  be  exploited  to  rebuild  or  ex¬ 
pand  an  RM  without  requiring  re-measurement.  New  RMs  can  be  constructed  from 
a  single  phase  map  with  comparable  and  even  improved  enhancement  performance. 
Since  this  process  is  linear,  a  similar  relationship  is  expected  for  the  focal-plane  sys¬ 
tem  using  the  inverse  spatial  Fourier  transform  of  the  tilt  phasemaps  in  Figure  16. 
However,  development  of  this  relationship  is  left  for  future  work.  This  capability  will 
be  needed  for  the  intended  application  of  imaging  around  corners.  Using  the  imaging 
system,  a  RM  can  be  constructed  and  then  extended  into  an  occluded,  unmeasured 
portion  of  the  scene.  The  scene  is  then  illuminated  using  reflective  inverse  diffusion 
via  the  RM  and  the  occluded  scene  can  be  imaged  from  the  point  of  view  of  the 
reflective  surface  using  dual  photography  [2], 
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7.1  Contributions 


A  collabortive  effort  with  the  previous  AFIT  masters  student,  Lieutenant  Jessica 
Schafer  (now  Captain  Jessica  Ullom)  [38,42]  who  collected  the  original  experimental 
reflective  inverse  diffusion  data  was  combined  with  the  corroborative  mathemati¬ 
cal  foundation  and  simulation  results  presented  in  Chapter  III  to  produced  the  first 
peer  reviewed  paper  on  the  subject,  “Reflective  Inverse  Diffusion,”  by  Burgi,  Ullom, 
Marciniak  and  Oxley.  The  paper  has  been  accepted  for  publication  in  Applied  Sci¬ 
ences  [9]. 

The  initial  results  in  the  development  of  RMs  for  reflective  inverse  diffusion  dis¬ 
cussed  in  Chapter  IV  was  presented  at  the  Optics  and  Photonics  conference  in  San 
Diego,  California  in  August  2016.  The  corresponding  SPIE  proceedings  paper  “Ma¬ 
trix  methods  for  reflective  inverse  diffusion,”  by  Burgi,  Marciniak,  Nauyoks  and  Oxley 
was  published  in  SPIE  Proceedings  volume  9961  [43]. 

The  matrix  methods  were  expanded  to  include  changes  to  the  optical  system  and 
material  surface  properties  that  can  affect  enhancement  performance  of  the  RM  in 
Chapter  V.  The  paper  “Measuring  the  reflection  matrix  of  a  rough  surface,”  by 
Burgi,  Marciniak,  Nauyoks  and  Oxley  has  been  submitted  for  publication  in  Optical 
Engineering  and  is  currently  under  peer  review  [44], 

A  final  paper,  “Expansion  of  the  Reflection  Matrix  of  a  Rough  Surface,”  by  Burgi, 
Marciniak,  Nauyoks  and  Oxley,  based  on  the  work  in  Chapter  VI  is  currently  in 
preparation  and  targeted  for  submission  to  Optical  Engineering. 

7.2  Future  Research 

Further  research  is  required  to  solidify  the  RM  measurement  method.  Improving 
the  vibration  isolation  of  the  optical  system  will  be  key  to  improving  the  RM  enhance¬ 
ment  performance.  Mathematical  methods  for  combating  the  effects  of  vibration,  such 
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as  dynamic  data  driven  applications  systems  (DDDAS)  [40],  can  now  be  explored  as 
a  more  cost  effective  and  suitable  solution  based  on  the  intended  application. 

Development  of  an  application-specific  optical  system  should  be  pursued  in  the 
future  to  replace  the  simple  single  lens  providing  demagnification.  Enhancement  per¬ 
formance  is  expected  to  increase  substantially  with  improved  optics.  An  imaging 
system  with  greater  levels  of  demagnification  should  improve  enhancement  perfor¬ 
mance  for  rough  surfaces  with  smaller  lateral  surface  dimensions  (£w2). 

The  substantial  RMS  phase  difference  between  the  measured  RM  and  a  the  RM 
reconstructed  from  a  single  phase  map,  suggests  multiple  RMs  exist  for  a  given  rough 
surface  with  comparable  enhancement  performances.  The  parallel  wavefront  modu¬ 
lation  method  used  for  measuring  the  RM  should  be  examined  more  closely  for  the 
cause  of  the  multiple  solutions.  Currently,  the  parallel  wavefront  modulation  opti¬ 
mizes  the  intensity  of  a  single  observation-plane  segment.  Future  research  should  also 
look  for  ways  to  optimize  the  RM  measurement  method  to  produce  the  maximum 
enhancement. 

Also,  currently,  measurement  of  the  RM  requires  less  time  and  provides  more 
information  that  iterative  optimization  techniques.  However,  with  the  linear  phase 
relationship  of  the  imaging  system  the  RMs  can  be  constructed  around  phase  maps  op¬ 
timized  with  more  sophisticated  intensity-based  phase-retrieval  methods  (e.g.,  steep¬ 
est  descent,  Gerchberg-Saxton,  conjugate-gradient  method)  and  should  be  examined 
for  practical  implementation  in  future  research.  The  RM  reconstruction  methods  for 
the  focal-plane  system  based  on  the  inverse  2-D  spatial  Fourier  transform  of  the  linear 
tilt  phase  map  should  be  also  be  investigated  in  future  work  to  determine  its  potential 
for  use  in  the  imaging-around-corners  application. 

The  coherence  length  of  the  laser  source  may  become  a  concern  at  the  longer 
propagation  distances  given  the  intended  application.  The  temporal  coherence  of  the 
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light  from  each  of  the  SLM  segments  decreases  as  the  propagation  distance  exceeds  the 
coherence  length  of  the  source.  As  the  degree  of  coherence  decreases,  it  will  become 
increasingly  difficult  for  the  SLM  to  compensate  for  the  surface  height  fluctuations  of 
the  reflector  and  produce  a  converging  phase  front  after  reflection.  Currently  available 
stabilized  HeNe  lasers  have  coherence  lengths  of  several  hundred  meters. 

There  are  significant  device  limitations  to  be  overcome  before  reflective  inverse 
diffusion  can  be  incorporated  into  a  system  for  imaging  around  corners.  High-speed 
phase  modulators  are  needed  for  the  intended  application.  Dual  photography  requires 
1024  intensity  measurements  to  produce  a  single  low-resolution  32  x  32  pixel  image 
from  reflected  light  [2].  Imaging  a  dynamic  scene  at  just  1  frame-per-second  would 
require  a  phase  modulator  with  a  frequency  response  of  approximately  1  kHz,  well 
beyond  the  the  30  Hz  of  the  current  LCoS  SLMs.  High  speed  phase  modulators 
would  also  reduce  the  measurement  time  of  the  RM  and  make  it  less  susceptible 
to  the  vibrational  noise  responsible  lower  than  expected  enhancement  performance. 
However,  these  speeds  would  likely  have  to  exceed  10’s  of  kHz. 
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Appendix  A.  MATLAB®  Source  Code 

A.l  Reflective  Inverse  Diffusion  -  Iterative  Algorithm 

The  following  MATLAB®  code  was  use  to  simulate  the  original  proof-of-concept 
experiments  from  Chapter  II  using  the  diffraction  base  model  from  Chapter  III.  The 
source  code  used  to  propagate  the  SLM  to  the  observation  plane  is  available  in  section 
A. 3.1. 


%  number  for  physical  SLM  pixels 
%  zero-pad  factor 


%%  Propagation  Simulations 
lambda=633*10'-9;  %  HeNe  wavelength 

k=2*pi/lambda;  %  wavenumber 

slmpix=512 ; 
a=2 ; 

D=0 . 00768; 
dxl=D/slmpix ; 
xl= ( (-slmpix/ 2 )  :  (slmpix/ 2-1 ) ) *dxl ;  % 

flens=0.5;  % 

zl=0 . 15 ;  % 

z2=0 .5;  % 

dx2=lambda*f lens/ (a*D) ; 

x2= ( (-a* slmpix/ 2 ) : (a* slmpix/ 2-1 ) ) *dx2; 


physcial  SLM  side  length 

SLM  delta-x 

SLM  coordinate  array 

positive  lens  focal  length 

distance  from  the  SLM  to  the  lens 

distance  from  the  reflector  to  the  CCD 

%  sample  pixel  size 
%  sample  coordinate  array 


sample=exp (li*2*pi*rand (a*slmpix) ) ;  %  uniform  distribution 

CCDchannels=1024;  %  total  CCD  segments 

CCDdimension=sqrt (CCDchannels)  ; 

SLMchannels=1024 ;  %  final  number  of  optimized  SLM  segments 

SLMdimension=sqrt ( SLMchannels ) ; 


sl=l; s2=512; s3=sl; s4=s2; 


%%  Iterative  Optimization  to  TestChannel 
%  start  with  random  SLM  segments 

SLMl=ExpImage (exp (li*2*pi*rand (slmpix) ) , [slmpix  slmpix]); 

TestChannel= (CCDdimension+1 ) *CCDdimension/2 ;  %  select  target  channel 

[ChannelMask, BackgroundMask] =ChannelMasks (TestChannel,  CCDdimension,  (s2-sl  +  l) )  ; 
PhaseSpacing=l 6;  %  number  of  tested  phase  values 

PhaseValue=exp (li*2*pi/PhaseSpacing* ( (-PhaseSpacing/2) : (PhaseSpacing/2-1) ) ) ; 
kmax=log ( SLMchannels ) /log ( 4 ) ;  %  calculate  iterations  of  outer  for  loop 

for  k=l:kmax 

N=2"(2*k);  %  number  of  segments  to  pre-optimize 

for  n=l : N 

ChannelIntensity=zeros (length (PhaseValue) , 1) ;  %  initializes  intensity 

SLMmask=ChannelMasks (n, sqrt (N) , slmpix) ;  %  select  SLM  segment 

for  p=l : length (PhaseValue) 
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tested  phase  value 


SLM1 (SLMmask) =PhaseValue (p) ; 

%  propagate  to  Obs  plane 

Uobs=PropagateSLM (SLM1,  zl,  flens,  z2,  sample,  a)  ; 

Iobs=abs (Uobs ( si : s2  ,  s3 : s4 ) )  . "2 ;  %  calculate  intensity 

%  save  intensity 

Channellntensity (p) =sqrt (mean (lobs (ChannelMask)  . ~2) )  ; 

end 

[~, idx] =max (Channellntensity )  ;  %  determine  maximum  intensity 

SLM1 (SLMmask) =PhaseValue (idx) ;  %  set  segment  to  phase  value  of  max 

end 

end 


A. 2  Reflective  Inverse  Diffusion  -  Matrix  Method 
A. 2.1  Focal  Plane  Simulation. 

The  following  MATLAB®  code  was  use  to  simulate  the  measurement  of  the  re¬ 
flection  matrix  (RM)  using  the  optical  setup  from  Figure  2  with  the  rough  surface 
reflector  placed  at  the  lens  focus.  The  diffraction  base  model  from  Chapter  III  was 
used.  The  source  code  used  to  propagate  the  SLM  to  the  observation  plane  is  available 
in  section  A. 3.1. 


%  number  for  physical  SLM  pixels 
%  zero-pad  factor 


%%  Propagation  Simulations 
lambda=633*10'-9;  %  HeNe  wavelength 

k=2*pi/lambda;  %  wavenumber 

slmpix=512 ; 
a=2 ; 

D=0 . 00768; 
dxl=D/slmpix ; 
xl= ( (-slmpix/ 2 )  :  (slmpix/ 2-1 ) ) *dxl ;  % 

flens=0.5;  % 

z 1=0 . 15 ;  % 

z2=0.5;  % 

dx2=lambda*f lens/ (a*D) ; 

x2= ( (-a*slmpix/2 ) : (a*slmpix/2-l) ) *dx2; 


physcial  SLM  side  length 

SLM  delta-x 

SLM  coordinate  array 

positive  lens  focal  length 

distance  from  the  SLM  to  the  lens 

distance  from  the  reflector  to  the  CCD 

%  sample  pixel  size 
%  sample  coordinate  array 


%%  Random  Sample 

SampleDimension=2048;  %  sample  size  in  pixels 
JitterSpace=0 ;  %  magnitude  of  random  shift 

JitterOf fset=l;  %  offset  center  for  random  shift 

%  sample  with  uniform  phase  distribution 

Jsample=exp ( li*2*pi*rand (SampleDimension+2* JitterSpace) ) ; 


%%  Device  Error 

CCDerror=0;  %  intensity  measurement  error  in  percent 

SLMerror=0;  %  SLM  error  in  radian 
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%%  Parallel  Wavefront  Opt  Variables 
CCDchannels=128~2;  %  total  CCD  segments 
CCDdimension=sqrt (CCDchannels)  ; 

CCDmap=re shape ( 1 : CCDchannels , CCDdimension, CCDdimension) ; 

SLMchannels=1024 ;  %  total  SLM  segments 
SLMdimension=sqrt ( SLMchannels )  ; 

G=SLMchannels/ 2 ; 

P=1:G; 

w0=2*pi;  %  1  Hz  frame  per  second 

wp=(G+p)/ (4*G)*wO;  %  SLM  segment  frequencies 

mapl=flip ( (checkerboard (1, SLMdimension/2 , SLMdimension/2 ) >0 .5) ,2) ; 
grpl=zeros (SLMdimension) ; 

grpl (mapl ) =wp;  %  group  1  checkboard  of  modulated  and  static  segments 

map2= (checkerboard (1, SLMdimension/2, SLMdimension/2) >0 . 5) ; 
grp2=zeros (SLMdimension) ; 

grp2 (map2 ) =wp;  %  group  2  checkboard  of  modulated  and  static  segments 
%%  Observation  Plane  Windowing 

sl=129; s2=896; s3=sl; s4=s2;  %  for  standard  20.6  micro  768x768  obs  plane 

%  sl=257 ; s2=17 92 ; s3=sl ; s4=s2 ;  %  for  standard  10.3  1536x1536  obs  plane 

%%  Initialize  Reflection  Matrix 
Rm=zeros (CCDchannels, SLMchannels)  ; 

%%  capture  group  1 
frames=zeros (CCDchannels,  4*G)  ; 
for  t  =  0 :  ( 4 *G-1 ) 

SLMnoise=SLMerror * (2*rand (size (grpl) ) -1) ; 

SLM=ExpImage (exp (li*quant (grpl*t+SLMnoise, 2*pi/16384) ) , [slmpix  slmpix] ) ; 

%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpace+1 ) - ( JitterSpace+1 ) ; 

Jy=randi (2* JitterSpace+1 ) - (JitterSpace+1) ; 

JstartX= Jitter Space! Jx+ Jit ter Of f set ; 

JstopX=JstartX+SampleDimension-l ; 

Jst art Y= Jitter Space! Jy+ Jit ter Of f set ; 

JstopY=JstartY!SampleDimens ion-1 ; 

sample=Js ample ( JstartX : JstopX,  JstartY : JstopY)  ; 

[Uobs, Uminus, Uplus] =PropagateSLM (SLM, zl, flens, z2, sample, a) ; 
CCDnoise=CCDerror* (2*rand (size(Uobs(sl:s2,sl:s2) ) ) -1) ll; 

Iobs=CCDnoise .*abs(Uobs(sl:s2,s3:s4)) .~2; 
frames ( : , ttl ) =re shape (Bin Image ( lobs , . . . 

[CCDdimension  CCDdimension] ) , CCDchannels, 1) ; 

end 

%%  Extract  Group  1  <ym|R>*tmp 
fftgrp=zeros (size (frames) )  ; 
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hmapl=re shape (mapl . ' , 1 , SLMchannels ) ; 
rmapl=logical (ones (CCDchannels , 1) *mapl 
fftgrp=fft (frames,  [  ] , 2) ; 

Rm (rmapl ) =f ft grp ( : , (G+2 ) : (2*Gll ) ) ; 

%%  capture  group  2  intensity  patterns 
frames=zeros (CCDchannels,  4*G)  ; 
for  t=0 :  ( 4 *G-1 ) 

SLMnoise=SLMerror * (2*rand (size (grp2) ) -1)  ; 

SLM=ExpImage (exp (li*quant (grp2*t!SLMnoise, 2*pi/16384) ) , [slmpix  slmpix] ) ; 

%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpace+1 ) - ( JitterSpace+1 )  ; 

Jy=randi (2* JitterSpace+1 ) - (JitterSpace+1)  ; 

JstartX= Jitter Space! Jx+ Jit ter Of f set ; 

JstopX=JstartX!SampleDimension-l ; 

Jst art Y= Jitter Space! Jy+ Jit ter Of f set ; 

JstopY=JstartY!SampleDimension-l ; 

sample=Js ample ( JstartX : JstopX,  JstartY : JstopY)  ; 

[Uobs , Uminus , Uplus ] =PropagateSLM (SLM, zl, flens, z2, sample, a) ; 
CCDnoise=CCDerror* (2*rand (size (Uobs (si : s2, si : s2) ) ) -1) ll; 

Iobs=CCDnoise .*abs(Uobs(sl:s2,s3:s4))  .  ~2; 
frames ( : , til ) =re shape (Bin Image ( lobs , . . . 

[CCDdimension  CCDdimension] ) , CCDchannels, 1) ; 

end 

%%  Extract  Group  2  <ym|R>*tmp 

fftgrp=zeros (size (frames) )  ; 
hmap2=reshape (map2 . ' , 1 , SLMchannels ) ; 
rmap2=logical (ones (CCDchannels, 1) *map2 (:).'); 
fftgrp=fft (frames,  [  ]  , 2) ; 

Rm (rmap2 ) =f ft grp ( : , (Gl2 ) : (2*Gf 1 ) ) ; 

%%  capture  group  3  intensity  patterns 
grp3=2*pi* (mapl* (3/8) fmap2* (4/8) ) ; 
frames=zeros (CCDchannels,  8)  ; 
for  t  =  0 : 7 

SLMnoise=SLMerror * (2*rand (size (grp3) ) -1) ; 

SLM=ExpImage (exp ( li*quant (grp3*t!SLMnoise,  2*pi/16384) )  ,  [slmpix  slmpix] )  ; 

%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpacell ) - ( JitterSpacef 1 ) ; 

Jy=randi (2* JitterSpacell ) - (JitterSpacell) ; 

JstartX= Jitter Space! Jxl JitterOf f set ; 

JstopX=JstartXlSampleDimension-l; 

Jst art Y= Jitter Space! Jyl JitterOf f set ; 

JstopY=JstartYl SampleDimension-1 ; 

sample=Js ample (JstartX : JstopX, JstartY : JstopY) ; 

[Uobs , Uminus , Uplus ] =PropagateSLM (SLM, zl,  flens,  z2,  sample,  a)  ; 
CCDnoise=CCDerror * (2*rand (size(Uobs(sl:s2,sl:s2) ) ) — 1) ll; 
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Iobs=CCDnoise .*abs(Uobs(sl:s2,s3:s4)) ."2; 
frames ( : , t+1 ) =re shape (Bin Image ( lobs , . . . 

[CCDdimension  CCDdimension] ) , CCDchannels,  1)  ; 

end 

%%  Calculate  Compenstation  Term  from  Group  3 
fftgrp=zeros (size (frames) ) ; 
fftgrp=fft (frames,  [  ] , 2) ; 

Ctemp=f ft grp ( : , 4 ) . *f ft grp ( : , 5 ) . /f ft grp ( : , 8 ) +f ft grp ( : , 4 ) +f ft grp ( : , 5 ) +f ft grp ( : , 8 ) ; 
C=Ctemp . / abs (Ctemp) ; 

Cm=repmat (C,  1 ,  G)  ; 

%%  Apply  Compensation  term 
Rcomp=Rm; 

Rcomp (rmap2 ) =Rcomp ( rmap2 ) . *Cm ( : ) ; 


A. 2. 2  Image  Plane  Simulation. 

The  following  MATLAB®  code  was  used  to  simulate  the  measurement  of  the 
reflection  matrix  (RM)  using  the  optical  setup  from  Figure  10  with  the  rough  surface 
reflector  in  an  image  plane  of  the  SLM.  The  diffraction  base  model  from  Chapter 
III  was  used;  however,  equation  (52)  provided  the  output  for  the  observation  plane. 
The  source  code  used  to  propagate  the  SLM  to  the  observation  plane  is  available  in 
section  A. 3. 2. 


%%  Propagation  Simulations 
lambda=633*10,'-9;  %  HeNe  wavelength 

k=2*pi/lambda;  %  wavenumber 

zl=0.5;  %  distance  from  sample  to  CCD 

%%  Programmable  variables 

Demagnif ication=7 . 68 ;  %  ideal  demagnification  of  SLM 

SLMdimension=32 ;  %  SLM  segments  per  side 

%%  Random  Sample 

SampleDimension=2048;  %  sample  size  in  pixels 
JitterSpace=0 ;  %  magnitude  of  random  shift 

JitterOf fset=l;  %  offset  center  for  random  shift 

%  sample  with  uniform  phase  distribution 

Jsample=exp ( li*2*pi*rand (SampleDimension+2* JitterSpace) ) ; 


CCDerror=0 ; 
SLMerror=0 ; 


%  intensity  measurement  error  in  percent 
%  SLM  error  in  radian 


%%  Virtual  SLM 

D=0 . 007  68/Demagnification;  %  SLM  side  length; 


dxl=D/SLMdimension; 

dx2=D/SampleDimension; 


%  SLM  pixel  size  in  meters 
%  pixel  size  in  the  observation  plane 
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%  Parallel  Wavefront  Optimization 
SLMchannels=SLMdimension''2 ; 
G=SLMchannels/ 2 ; 
p=l:G; 
w0=2  *pi ; 

wp= (G+p) / (4*G) *w0; 


Parameters 

%  Total  SLM  segments 
%  SLM  segments  per  group 
%  array  counter 

%  1  Hz  frame  per  second  sampling 
%  SLM  segment  frequencies 


mapl=flip ( (checkerboard (1, SLMdimension/2 , SLMdimension/2 ) >0 .5) ,2) ; 
grpl=zeros (SLMdimension) ; 

grpl (mapl ) =wp;  %  group  1  checkboard  of  modulated  and  static  segments 


map2= (checkerboard (1, SLMdimension/2, SLMdimension/2) >0 . 5) ; 
grp2  =  zeros (SLMdimension)  ; 

grp2 (map2 ) =wp;  %  group  2  checkboard  of  modulated  and  static  segments 


%%  Observation  plane  windowing 

if  SampleDimension==32 

sl=25 ; s2=344 ;  %  for  sample  32x32 

ScanStep=l ; 

elseif  SampleDimension==64 

sl=369; s2=1008;  %  for  sample  64x64 

ScanStep=l ; 

elseif  SampleDimension==96 

sl=481; s2=2592;  %  for  sample  96x96 

ScanStep=l ; 

elseif  SampleDimension==128 

sl=2049; s2=3328;  %  for  sample  128x128 

ScanStep=l ; 

else 

error ( ' Unsupported  SampleDimension ' ) 

end 


CCDdimension=32 ;  %  obs  plane  width  in  segments 

CCDchannels=CCDdimension~2;  %  total  number  of  CCD  segments 
CCDmap=reshape (1 : CCDchannels, CCDdimension, CCDdimension) ; 

%%  Initialize  Reflection  Matrix 
Rm=zeros (CCDchannels, SLMchannels)  ; 

%%  capture  group  1 
frames=zeros (CCDchannels,  4*G)  ; 
for  t  =  0 :  ( 4 *G-1 ) 

SLMnoise=SLMerror* (2*rand (size (grpl) ) -1) ; 

SLM=exp ( li*quant (grpl*t+SLMnoise, 2*pi/16384) ) ; 

%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpace+1 ) - ( JitterSpace+1 )  ; 

Jy=randi (2* JitterSpace+1 ) - (JitterSpace+1)  ; 

JstartX= Jitter Space! Jx+ Jit ter Of f set ; 
JstopX=JstartX+SampleDimension-l; 

Jst art Y= Jitter Space! Jy+ Jit ter Of f set ; 
JstopY=JstartY!SampleDimension-l ; 
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sample=Js ample ( JstartX : JstopX, JstartY : JstopY) ; 


Uobs=PropSLMImage (SLM, D, sample, zl )  ; 

CCDnoise=CCDerror* (2*rand (size(Uobs(sl:s2,sl:s2) ) ) -1) +1; 
Iobs=CCDnoise .*abs(Uobs(sl:s2,sl:s2)) .  ~2; 
frames ( : , t  +  1 ) =re shape (Bin Image ( lobs ,  .  .  . 

[CCDdimension  CCDdimension] ) , CCDchannels, 1) ; 

end 

%%  Extract  Group  1  <ym|R>*tmp 
fftgrp=zeros (size (frames) )  ; 
hmapl=re shape (mapl .  '  ,  1 ,  SLMchannels )  ; 
rmapl=logical (ones (CCDchannels, 1) *mapl (:).'); 
fftgrp=fft (frames,  [  ] , 2) ; 

Rm (rmapl ) =f ft grp ( : , (G+2 ) : (2*G+1 ) ) ; 

%%  capture  group  2  intensity  patterns 
frames=zeros (CCDchannels,  4*G)  ; 
for  t  =  0 :  ( 4 *G-1 ) 

SLMnoise=SLMerror * (2*rand (size (grp2) ) -1) ; 

SLM=exp ( li*quant (grp2*t+SLMnoise, 2*pi/16384) ) ; 

%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpace+1 ) - ( JitterSpace+1 )  ; 

Jy=randi (2* JitterSpace+1 ) - (JitterSpace+1)  ; 
JstartX=JitterSpace+ Jx+ Jit ter Of f set ; 
JstopX=JstartX+SampleDimension-l ; 

JstartY=JitterSpace+ Jy+ Jit ter Of f set ; 

Jst opY=Jst art Y+SampleDimens ion-1 ; 

sample=Js ample (JstartX : JstopX, JstartY : JstopY) ; 

Uobs=PropSLMImage (SLM, D, sample, zl) ; 

CCDnoise=CCDerror* (2*rand (size(Uobs(sl:s2,sl:s2) ) ) -1) +1; 
Iobs=CCDnoise .*abs(Uobs(sl:s2,sl:s2)) . ~2 ; 
frames ( : , t+1 ) =re shape (Bin Image ( lobs , . . . 

[CCDdimension  CCDdimension] ) , CCDchannels, 1) ; 

end 

%%  Extract  Group  2  <ym|R>*tmp 

fftgrp=zeros (size (frames) )  ; 
hmap2=reshape (map2 . ' , 1 , SLMchannels ) ; 
rmap2=logical (ones (CCDchannels, 1) *map2 (:).'); 
fftgrp=fft (frames,  [  ]  , 2) ; 

Rm (rmap2 ) =f ft grp ( : , (G+2 ) : (2*G+1 ) ) ; 

%%  capture  group  3  intensity  patterns 
grp3=2*pi* (mapl* (3/8) +map2* (4/8) ) ; 
frames=zeros (CCDchannels,  8)  ; 
for  t  =  0 : 7 

SLMnoise=SLMerror * (2*rand (size (grp3) ) -1) ; 

SLM=exp ( li*quant (grp3*t+SLMnoise, 2*pi/16384) ) ; 
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%  simulation  optics  table  vibration 
Jx=randi (2* JitterSpace+1 ) - ( JitterSpace+1 )  ; 

Jy=randi (2* JitterSpace+1 ) - (JitterSpace+1) ; 

JstartX= Jitter Space! Jx+ Jit ter Of f set ; 

JstopX=JstartX+SampleDimension-l ; 

Jst art Y= Jitter Space! Jy+ Jit ter Of f set ; 

Jst opY=Jst art Y!SampleDimens ion-1; 

sample=Js ample ( JstartX : JstopX, JstartY : JstopY) ; 

Uobs=PropSLMImage (SLM, D, sample,  zl )  ; 

CCDnoise=CCDerror* (2*rand (size(Uobs(sl:s2,sl:s2) ) ) -1) !l; 

Iobs=CCDnoise .*abs(Uobs(sl:s2,sl:s2)) .  ~2; 
frames ( : , til ) =re shape (Bin Image ( lobs ,  .  .  . 

[CCDdimension  CCDdimension] ) , CCDchannels,  1)  ; 

end 

%%  Calculate  Compenstation  Term  from  Group  3 
fftgrp=zeros (size (frames) )  ; 
fftgrp=fft (frames,  [  ] , 2) ; 

Ctemp=f ft grp ( : , 4 ) . *f ft grp ( : , 5 ) . /f ft grp ( : , 8 ) if ft grp ( : , 4 ) if ft grp ( : , 5 ) if ft grp ( : , 8 ) ; 
C=Ctemp . / abs (Ctemp) ; 

Cm=repmat (C,  1 ,  G)  ; 

%%  Apply  Compensation  term 
Rcomp=Rm; 

Rcomp (rmap2 ) =Rcomp ( rmap2 ) . *Cm ( : ) ; 


A. 3  SLM  Propagation  Functions 
A. 3.1  Focal  Plane  Optical  System. 

The  following  MATLAB®  code  was  used  to  propagate  the  simulated  SLM  to  the 
observation  plane  base  on  the  focal  plane  optical  setup  in  Figure  2  with  the  rough 
surface  reflector  placed  at  the  lens  focus.  The  diffraction  base  model  from  III  was 
used.  The  propagation  was  based  on  equation  (38). 

%  Uobs  =  observation  plane 

%  Uminus  =  field  just  prior  to  reflection 
%  Uplux  =  field  just  after  reflection 
%  dx3  =  observation  plane  pixel  size 
%  zl  =  distance  from  SLM  to  Lens 
%  z2  =  distance  from  Lens  to  Observation  plane 

function  [Uobs , Uminus , Uplus , dx3 ]  = 

Propagate SLM ( input , zl, focus, z2, sample,  a,  re calc) 

%  persisent  variable  for  computational  efficiency 
persistent  prevZl  prevF  prevZ2  prevM  prevN  prevA 
persistent  x3  chirp  dxl 
persistent  Hslm21ens  Hsample2ccd 
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if  nargin<7 
recalc=0 ; 
elseif  nargin<6 
a=2  ; 

recalc=0 ; 
elseif  nargin<5 
sample=l ; 
a=2 ; 

recalc=0 ; 

end 

[m, n] =size (input) ; 

M=a*m; 

[msample, nsample] =size (sample) ; 
if  msample~=l  | |  nsample~=l 

if  a*m~=msample  &&  a*n~=nsample 

error (' sample  incorrect  size') 

end 

end 

%  check  for  changes  that  require  recalculation  of  persistent  variables 
if  isempty (prevZl) 
recalc=l ; 
prevZl=0 ; 
prevF=0 ; 
prevZ2=0 ; 
prevM=0 ; 
prevN=0 ; 
prevA=0 ; 

elseif  zl~=prevZl  |  |  focus~=prevF  |  |  ,  .  . 

z2~=prevZ2  |  |  m~=prevM  | |  n~=prevN  |  |  a~=prevA 
recalc=l ; 

end 

if  recalc 

display (' recalculating '  )  ; 

lambda=633*10  ~-9; 
k=2 *pi/ lambda  ; 

D=0. 00768;  %  SLM  dimension 

Ll=a*D; 

%  slm  plane  coordinates 

dxl=Ll/M;  %  SLM  pixel  size 

xl=-Ll/2 : dxl : Ll/2-dxl ; 

[XI , Y1 ] =meshgrid (xl )  ; 

%  propagation  TF  -  SLM  to  Lens 
fx2= (-M/2 : M/2-1) *1/L1; 


91 


[FX2 , FY2 ] =meshgrid ( f x2 )  ; 

Hslm21ens=exp (li*k*zl*sqrt (1- (lambda*FX2) . ~2- (lambda*FY2) . "2) ) ; 

%  sample  plane  coordinates 
dx3=lambda*focus/Ll ; 

L3=lambda*focus/ dxl ; 
x3=-L3/2 :dx3 :L3/ 2-dx3; 

[X3 , Y3 ] =meshgrid (x3 ) ; 

%  Fraunhofer  chirp 

chirp=exp (li*k/ (2*focus) * (X3. ~2+Y3 . ~2) ) ; 

%  propagation  TF  -  sample  to  CCD 
f x4= (-M/2 :M/2-l) *1/1,3; 

[FX4  ,  FY4 ] =meshgrid ( f x4 ) ; 

Hsample2ccd=exp ( li*k*z2*sqrt (1- (lambda*FX4) . "2- (lambda*FY4) . ~ 2) ) ; 

prevZl=zl ; 
prevF=f ocus ; 
prevZ2=z2 ; 
prevM=m; 
prevN=n ; 
prevA=a; 

end 

ul=padarray ( input , [ (M-m) 12,  (M-m) / 2 ] ) ; 

Uminus=chirp  ,*fftshift  (fft2  (fftshift  (ul)  )  )  *dxl,'2.  *Hslm21ens ; 
Uplus=sample . *Uminus ; 

Uobs=f ftshift (ifft2 (fft2 (fftshift (Uplus) ) . *fftshift (Hsample2ccd) ) ) ; 
end 


A. 3. 2  Image  Plane  Optical  System. 

The  following  MATLAB®  code  was  used  to  propagate  the  simulated  SLM  to  the 
observation  plane  base  on  the  imaging  system  optical  in  Figure  10  with  the  rough 
surface  reflector  placed  is  the  imaging  plane  of  the  SLM.  The  diffraction  base  model 
from  III  was  used.  The  propagation  was  based  on  equation  (52). 


%  ObsOut  =  observation  plane 

%  zl  =  distance  from  sample  to  Observation  plane 

function  [ObsOut]  =  PropSLMImage (SLM1, SLMdim, sample, PropZ, recalc) 

%  persistent  variables  for  computational  efficiency 
persistent  pad  pSLMdim  pZ  Hsample2ccd 
persistent  pSx  pSLMx 

[SLMx, SLMy ] =size (SLM1) ; 

[Sx, Sy]=size (sample)  ; 

if  nargin<5 


92 


recalc=0;  %  default  use  previous  Hsample2ccd  values 

end 

if  isempty (Hsample2ccd) 
recalc=l ; 

elseif  SLMdim~=pSLMdim  | |  PropZ~=pZ  | |  Sx~=pSx  | |  SLMx~=pSLMx 
recalc=l ; 

end 

if  mod ( Sx, SLMx) ~=0  ||  mod (Sy, SLMy) ~=0 

error (' sample  size  not  multiple  of  SLM') 

elseif  Sx~=Sy  | |  SLMx~=SLMy 

error (' sample  or  SLM  not  square') 

elseif  mod (SLMx, 2 ) ~=0  ||  mod(Sx,2)~=0 

error ('even  samples  only') 

end 

if  recalc 

lambda=633*10  ~-9; 
k=2 *pi/ lambda; 


if  Sx==32 

pixout=3  68 ; 
elseif  Sx==64 

pixout  =  137  6 ; 
elseif  Sx==96 

pixout=3072 ; 
elseif  Sx==128 
pixout=537  6 ; 

else 

error (' Unsupported  Sample  Size'); 

end 

%  side  length  of  demagnified  zero  padded  SLM  image 
Ll=pixout* (SLMdim/Sx) ; 

f x= (-pixout /2 : pixout / 2-1) *1/L1; 

[FX, FY] =meshgrid (fx) ; 

%  Rayleigh-Sommerfeld  Transfer  Function 

Hsample2ccd=exp ( li*k*PropZ*sqrt (1- (lambda*FX)  . ~2- (lambda*FY)  . ~ 2) )  ; 
%  Fresnel  Transfer  Function 

%Hsample2ccd=exp (-li*pi*lambda*PropZ* (FX. "2+FY. ~2) )  ; 

pad= (pixout-Sx) /2 ;  %  pad  output 

pSLMdim=SLMdim; 

pZ=PropZ ; 

pSx=Sx; 

pSLMx=SLMx; 

end 

ul=padarray ( sample . *ExpImage ( SLM1 ,[ Sx  Sy] ) , [pad  pad]); 

ObsOut=f ftshift (ifft2 (fft2 (fftshift (ul) ) . *fftshift (Hsample2ccd) ) ) ; 
end 
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